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PREFACE 


The effort required to improve the NASA/MSFC Global Reference Atmospheric Model 
(GRAM) was sponsored by the Earth Science and Applications Division, Space Science Labora- 
tory, NASA Marshall Space Flight Center, under the technical monitorship of Mr. Dale L. 
Johnson, the NASA COR. The modeling work was accomplished by the Georgia Institute of 
Technology, School of Earth and Atmospheric Sciences, under NASA grant NAG8-078. The effort 
was funded largely by NASA Headquarters Code M. New Technology, Inc. (NTI) coordinated the 
transfer of program, and program checkout, for use on the MSFC IBM computer. 


This report presents the various GRAM-90 program code sequences as configured to operate 
on MSFC’s IBM 3090 computer. Users should note that some of the code listings presented herein 
will require certain site-specific modifications to operate on other machines or configurations. 


The first three listings start with the GRAM-90 program FORTRAN code, followed by an 
example of the IBM Job Control Language (JCL) code required to run the program and an example 
of resulting output. The next two listings start with a combined JCL and FORTRAN code sequence 
used to convert an individual monthly four-dimensional (4-D) data set to the direct access format 
required by the GRAM-90, followed by an example of the values for the variables found in each 
data record. Since even | month of the 4-D data base is extremely large, any complete listing 
becomes infeasible, and thus is omitted from this volume. However, the last listing is an output of 
the entire SCIDAT9 data base, which, though much more extensive than its predecessor, is hereby 
included in accordance with previous complete GRAM updates. 


Key individuals who have either contributed toward the completion and testing of the 
GRAM-90 model (and/or report) are Steve Smith (MSFC) and Rhonda Blocker (BCSS). The tech- 
nical work and user help information has been documented in part I of this report. Qualified 
requestors may obtain (purchase) copies of the computer program for this NASA/MSFC GRAM by 
contacting: COSMIC, the University of Georgia, 382 E. Broad St., Athens, GA 30602, 
telephone (404) 542-3265. 
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TECHNICAL MEMORANDUM 


THE NASA/MSFC GLOBAL REFERENCE ATMOSPHERIC 
MODEL—1990 VERSION (GRAM-90) 


Part ll: Program/Data Listings 


1. GRAM-90 CODE 


The following is a listing of the 1990 Global Reference Atmosphere Model (GRAM-90) 
FORTRAN program code. As in previous versions of the GRAM code (Justus et al. 1980), each 
segment of the program (main program, subroutine, function) has lines which are uniquely identi- 
fied by a three- or four-character segment code coupled with a sequence number appearing at the 
extreme right of each line. 


PROGRAM GRAMSO GRAM 1 

C FIRST DATA LINE READS INITIAL HEIGHT (KM), INITIAL LATITUDE (OEG) GRAM 2 
Cc INITIAL LONGITUDE (DEG), F10.7. MEAN F10.7, AP, MONTH, DAY, GRAM 3 
C YEAR (TOTAL YEAR - 1900), GREENWICH HOUR, MINUTES, SECONDS, GRAM 4 
Cc LATITUDE INCREMENT (DEG), LONGITUDE INCREMENT (DEG), GRAM 5 
C HEIGHT DECREASE (KM), MAXIMUM NUMBER OF POSITIONS (EXCLUDING GRAM 6 
Cc INITIAL POSITION) TO BE COMPUTED, TIME INCREMENT BETWEEN GRAM 7 
Cc POSITIONS, TRAJECTORY OPTION, OUTPUT OPTION,MINIMUM GEOSTROPHICGRAM 8 
Cc LATITUDE GRAM 9 
COMMON /IOTEMP/IOTEM1, IOTEM2, [US,0D,xXMJD,PHIt,PHI, GRAM 10 

: NSAME .RP1, RD1, RT1, SP1, SOD1, ST1, RUT, RVI, SUT, SV1,GRAM 11 

%$ MN, IOA, TYR, Hi, PHI1R, THETIR,G,RI.H,PHIR, THETR,F10,F 108, AP, GRAM 12 

, IHR .MIN,NMORE ,DX.HL,VL.0Z,8,EPS, 1OPP,LOOK,1ET,.GLAT, GRAM 13 

IRP IS. RDIS.RTIS,RUIS,RV1S,.SP15,SD1S,ST1S,SU1S,SVIS, GRAM 14 
2U0S1,V0S1,U0L1,VDL1,U0$2, V0S2,UDL2,VDL2 GRAM 15 
COMMON /CHIC/LA(4,4),NB(2), IWSYM,USH, VSH,OUSH,OVSH GRAM 16 
COMMON /WINCOM/DUMSTF(14),UPRE,VPRE ,DUPRE ,DVPRE GRAM 17 

9090 FORMAT(‘1 ****** GLORAL REFERENCE ATMOSPHERIC MODEL - 1990’, GRAM 18 
& ' (GRAM9Q) *hrttt’/) GRAM 19 
PI=3.1415927 GRAM 20 
FAC=0.017453293 GRAM 21 
LOOK=0 GRAM 22 

MONTH = O GRAM 23 
IOPT=O GRAM 24 

H=0. GRAM 25 

5 IF(IOPT.£Q.0.0R.(IOPT.GT.O.AND.H.LT.O.)) GO TO 6 GRAM 26 
READ(IOPT,*})IET,H.PHI, THET GRAM 27 

GO TO 5 GRAM 28 

6 MN = MONTH GRAM 29 
NSAME = O GRAM 30 
READ(5.*,. END=90) H1.PHE!, THET1,F10,f 10B, AP MN, IDA, EYR,THRO,MINO, GRAM 31 

{ ISECO,OPHI ,.DTHET ,.DH,NMAX.INCT,IOPT,TOPP.GLAT GRAM 32 
WRITE(6,9090) GRAM 33 
IF(ABS(PHI1).LE.90.)GO TO 7 GRAM 234 

PHI 1=SIGN( 180. -ABS({PHI1),PHI1) GRAM 435 
THET4=THET141+180. GRAM 36 

IF (THET1.GT.360. )THET1=THET1-9360. GRAM 37 

7 IFC(THET1.0T.0.) THET1=THET1+*960 GRAM 38 
GLAT = ABS(GLAT) GRAM 39 

IF (GLAT.LT. 5S.) GLAT = 5. GRAM 40 

IF (GLAT.GE.18.) GLAT = 17.999 GRAM 41 
GLATF=GLAT*FAC GRAM 42 
WRITE(6.9010) H1,PHI1, THET1,F10,F10B,AP.MN.IDA ,IYR,ITHRO,.MINO, GRAM 43 

$ ISECO.OPHI ,OTHET,DH,NMAX,INCT, IOPT,FOPP,GLAT GRAM 44 

Cc SET NSAME TO AVOID SETUP GRAM 45 
15 IF (MN.EQ.MONTH) NSAME = 1 GRAM 46 

Cc LOOKUP ON MULTIPLE PASSES GRAM 47 
MONTH = MN GRAM 48 

Cc CONVERT LATITUDE TO RADIANS GRAM 49 
PHI 1R=PHI1*FAC GRAM 50 

C CONVERT LONGITUDE TO RADIANS GRAM 51 
THETIR=THET1*FAC GRAM 52 

Cc CONVERT LATITUDE INCREMENT TO RADIANS GRAM 53 


DPHIR=DPHI *F AC GRAM 54 


16 


17 


CONVERT LONGITUDE INCREMENT TO RADIANS 

DTHETR=DTHET*FAC 

READ DATA FILE TO INITIALIZE ARRAYS 

CALL SETUP 

NT = 4 

IF(IOPT.EQ.0) GO TO 18 

READ(IOPT,*)IET,H,PHI, THET 

IF(ABS(PHI).LE.90.) GO TO 16 

PHI = SIGN( 180. -ABS(PHI),PHI) 

THET = THET + 180. 

IF (THET.LT.O. )THET=THET+9360. 

IF(THET.GT.360.) THET = THET - 360. 

PHIR=PHI*FAC 

THE TR=THET*FAC 

GO tO 19 

H = H1 - OH 

DISPLACES POSITION BEFORE EVALUATION OF ATMOSPHERIC PARAMETERS 

IET = INCT 

PHIR=PHI 1R+DOPHIR 

THE TR=THETIR+DTHETR 

IF (ABS(PHIR).LE.PI/2.}) GO TO 17 

PHIR = SIGN(PI-ABS(PHIR),PHIR) 

DPHIR=-DPHIR 

THETR = THETR + PI 

IF (THETR.GT.2.*PI) THETR = THETR ~- 2.*PI 

IF (THETR.LT.O.) THETR = THETR + 2.*PI 

A=EQUATORIAL EARTH RADIUS, B=POLAR FARTH RADIUS 

EPS= EARTH ECCENTRICITY 

A = 6378.160 

B = 6356.7747 

EPS=(1.-(B*B)/(A*A)) 

COMPUTES RADIUS TO HEIGHT H, AND GRAVITY AT HEIGHT AND 
LATITUDE PHIR 

CALL RIG 

ISEC=ISECO+IET 

ISEC=MOD(ISEC,60) 


MIN = MINO + IET/60 
IHR = IHRO + MIN / 60 
MIN = MOO(MIN,60) 


COMPUTES P,0,T.U,V AT FIRST POSITION AFTER INITIAL POSITION 
IF(H1.LE.30.) LOOK=1 
IFC(ABS(PHIR).GT.GLATF) GO TO 195 
IF((H.GE.30.0).AND.(H.LE.85.0))GOTO 195 
PHI1S = PHI1R 
PHIS=PHIR 
DPHIS=(PHIR+GLATF)/(2.*GLATF ) 
PHIR=GLATF 
PHIIR = PHIR + PHI1S - PHIS 
CALL SCIMOD(0O) 
UP2=UPRE 
VP2=VPRE 
DUP 2=DUPRE 
DVP2=DVPRE 
PHIR=-GLATF 


GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 
GRAM 


CALL SCIMOD(O) GRAM 109 


UP 1=UPRE GRAM 110 
VP 1=VPRE GRAM 111 

DUP 1=DUPRE GRAM 112 

DVP 1=DVPRE GRAM 113 
UPRE=UP 1+(UP2-UP1)*DPHIS GRAM 114 
VPRE=VP14+(VP2-VP1)*DPHIS GRAM 115 

DUPRE =DUP 14+ (DUP2-DUP1)*DPHIS GRAM 116 
DVPRE=DVP 1+(DVP2-DOVP1)*DPHIS GRAM 117 

PHIIR = PHIIS GRAM 118 
PHIR=PHIS GRAM 119 

195 CALL SCIMOD(1) GRAM 120 
90 NT = NT + 1 GRAM 121 

IF (I1OPT.EQ.0) GO TO 22 GRAM 122 
READ(TOPT, *)IET,H,PHI, THET GRAM 123 
IF(H.LT.O.)GO TO 5 GRAM 124 
IF(ABS(PHI).LE.90.)GO TO 21 GRAM 125 
PHI=SIGN( 180. -ABS(PHI),PHI) GRAM 126 

THE T=THET+ #80. GRAM 127 

24 IF (THET.LT.O. )THET=THET*+360. GRAM 128 
IF(THET.GT.360. )THET=THET 360. GRAM 129 
PHIR=PHI*FAC GRAM 130 
THETR=THET*FAC GRAM 131 

GO TO 25 GRAM 132 

Cc INCREMENT THE HEIGHT GRAM 133 
22 H = Ht - DH GRAM 134 
IF (H .LT. 0.0) GO TO 5 GRAM 135 

C INCREMENT THE LATITUDE GRAM 136 
PHIR=PHIR+DPHIR GRAM 137 

C INCREMENT THE LONGITUDE GRAM 198 
THE TR=THETRIDTHETR GRAM 139 
Cle READS NEW INPUT IF ABS(LAT) GTR 90 DEG GRAM 140 
IF (ABS(PHIR).LE.PI/2) GO TO 23 GRAM 141 
PHIR=SIGN(PI-ABS(PHIR),PHIR) GRAM 142 
DPHIR=-DPHIR GRAM 143 
THETR=THETRIPI GRAM 144 

23 1F (THETR.GT.2.*P1) THETR = THETR - 2. * PI GRAM 145 

IF (THETR.LT.O.) THETR = THETR + 2. * PI GRAM 146 

C INCREMENT THE TIME GRAM 147 
IET=IET+INCT GRAM 148 

2956 MIN=MINO+IET/60 GRAM 149 
ISEC=ISECO+IET GRAM 150 
LSEC=MOD(ISEC,50) GRAM 151 
IHR=LHRO+MIN/6O GRAM 152 
MIN=MOD(MIN,60) GRAM 153 

C COMPUTE RADIUS AND GRAVITY AT NEW POSITION GRAM 154 
CALL RIG GRAM 155 

C COMPUTE P,0,T,.U,V, AT NEW POSITION GRAM 156 
IF(H1 .GE. 30.0 .AND. H .LT. 30.Q0)LOOK = ! GRAM 157 
IF(ABS(PHI).GT.75. .AND.ABS(PHI1).LE.75. )LOOK = 4 GRAM 158 

IF (ABS(PHIR).GT.GLATF) GO TO 80 GRAM 159 
IF((H.GE.30.0).AND.(H.LE.85.0))GOTO 80 GRAM 160 

PHI1S = PHIIR GRAM 161 


PHIS=PHIR GRAM 162 


OCG 


DPHIS=(PHIR+GLATFE )/(2.*GLATF ) GRAM 163 


PHIR=GLATF GRAM 164 
PHI1R = PHIR + PHI15 - PHIS GRAM 165 
CALL SCIMOD(O) GRAM 166 
UP2=UPRE | GRAM 167 
VP2=VPRE GRAM 168 
DUP 2=DUPRE GRAM 169 
DVP2=DVPRE GRAM 709 
PHIR=-GLATF | GRAM 110 
CALL SCIMOD(O) GRAM 172 
UP 1=UPRE GRAM 173 
VP 1=VPRE GRAM 174 
DUP 1=OUPRE GRAM 175 
DVP 1=OVPRE GRAM 176 
UPRE =UP 1+(UP2-UP1)*DPHIS | GRAM 177 
VPRE=VP1+(VP2-VP1)*DPHIS GRAM 178 
DUPRE =DUP 1+(DUP2-DUP1)*DPHIS GRAM 179 
DVPRE =DVP 1+(DOVP2-DVP1)*DPHIS GRAM 180 
PHIR=PHIS GRAM 181 
PHI1R = PHI1S GRAM 182 
80 CALL SCIMOD(1) GRAM 183 
GRAM 184 

Rio eee READS NEW INPUT IF NMORE = O OR MAX POINTS COMPUTED GRAM 185 
IF(NMORE .EQ.0.0R.(IOPT.EQ.0.AND.NT.GE.NMAX)) GO TO 5 GRAM 186 
CYCLE TO NEW POSITION GRAM 187 
GO TO 20 GRAM 188 
90 STOP GRAM 189 
9010 FORMAT(’ INITIAL HEIGHT = ’,F7.2,’ KM’,143,’‘INITIAL LAT = ’, GRAM 190 
1 FG.2,’ DEG’,1T83, INITIAL WEST LON = ’,F6.2,’ DEG’,/,’ F10.7 =’, GRAM 191 
2 F8.2,743,’MEAN F10.7 = ',F7.2,7T83,'AP = ‘,F8.2,/,’ DATE = ’, GRAM 192 
3 12,’/',12,‘/’,12,T43, ‘GREENWICH TIME = ’,12,':',12,’:’, 12/, GRAM 193 
4 ‘ LAT INCREMENT = ’,F6.2,’ DEG’ ,743,’WEST LON INCREMENT = °, GRAM 194 
5 F6.2,’ DEG’,T83, ’HEIGHT INCREMENT = ’,F7.2,’ KM’,/, GRAM 195 
G6 ‘ MAXIMUM NUMBER OF POSITIONS = ’,14,1T43,’TIME INCREMENT = °, GRAM 196 
7 14,’ SEC’ ,/2X,’* TRAJECTORY OPTION = ’,14,143, ‘OUTPUT OPTION = ’. GRAM 197 
8 12,783, ’MIN GEOSTROPH LAT = ’,F5.1,/) GRAM 198 
END GRAM 199 
COMMUN/C4/0DUN1(32),NG,P(16,26),0(16,26).7T(16,26),SP(16. 26) ADJU 2 
$,50(16,26),.$1(16,26),DU1,0U2.HS ADJU 3 
COMMON/ADUCOM/A(26,3), B(26), X(26). KOUNT ADJU 4 
DIMENSION PQ(26), QQ(26), UC(26), VC(26), wco(26), Uf26), V(26)., ADJU 5 
$ w(26) ADJU 6 
ASSUMPTIONS: ADJU 7 
HS IS THE SURFACE LEVEL ADJU B 
ALL DATA VALUES ABOVE SURFACE LEVEL ARE IN 1 KM INCREMENTS ADJU 9 
E1=0.075 ADJU 10 
E2=0.150 ADJUU 11 
MAXIT=3 ADJU 12 
KSMAX=10 ADJU 13 
HSU = HS ADJU 14 
IF (HS.LT.O.) HSJ = O. ADJU 15 
JJ=INT(HSJt2. ) | | ADJU 16 
1SS=1 ADJU 17 
CONST=28703./980.665 ADJU 18 
N=26 ADJU 19 


ITER=0 ADJU 20 


15 


12 


“fr = # © 


_ = © 2 © 


25 


90 B(I2)=UCIP)-U(I)-(U(T)-V(T)+Ww(T)) 


35 


2 «© # @ 


UC(1)=SORT(ABS(SP(KOUNT. 
VC(1)=SQRT(ABS(SD(KOUNT, 
wo(t)=SORT{ABS(ST(KOUNT, 
DO S I=dJ,N 
UC(1)=SORT(ABS(SP(KOUNT.1))) 
VC(1)=SORT(ABS(SD(KOUNT,1I)})) 
wo(1)=SQRT(ABS(ST(KOUNT,1))) 
NM=N- { 

NP=N+ 4 

SETS UP QUADRATURE FACTORS 
PQ(1)=500. *(FLOATCINT(HSJ+1. 


1))) 
1))) 
1))) 


))}-HS)/(CONST*T(KOUNT, 1)) 


Q0(1)=500. *(FLOAT(INT(HSu+ 1.) )-HS)/(CONST*T(KOUNT , JJ) ) 


DO 15 I=JJ.NM 

IP=I+1 
PQ(1)=500./(CONST*T(KOUNT,I1)) 
QQ(1)=500./(CONST*T(KOUNT,IP)) 
GO TO 58 

NM=N- t 

NP =N+ 4 

pO 14 1=1,26 

uU( I )=UCc(1)*UC(T) 
V(I)=vc(1)*vc (1) 
w(I}=we(1)*wc(1) 

CONTINUE 

INITIALIZE A(I.U) 

DO 20 1=1,26 

DO 20 J=1,3 

A(1,J)=0. 

SETS UP COEFFICIENTS 

12=0 

pO 35 I=1,NM 
IF(1.GT.t.AND.I.LT. JJ) GO TO 35 
Aw=1./SP(KOUNT.I) 
Bw=1./SD(KOUNT,I) 
Cw=1./ST(KOUNT,1) 


IM=I-1 
IF(I.EQ.JU) IMet 
IP=1+1 

IF (£.EQ.1) IP=dJd 
12=12+1 


AWt=1./SP(KOUNT,IP) 
Bwi=1./SO(KOUNT,IP) 
Cwt=1./ST(KOUNT,IP) 
IF(I1.EQ.1) GO TO 25 


A( 12. 1)=-(1.-QQ(IM))*(1.+PQ(1))/AWs (1. /BW+1./CW)*PQ(T) FQQ(TM) 
A(12.2)=((1.-QO(T))**2)/AW14( (1. +PQ(T))**2)/AW+( 1. /BWHt . /CW) 
$ *(PQ(1)**2)4(1. /BW1+1./CW1)}*OQ(1)**2 


IF(I.EQ.NM) GO TO 30 


A(12,3)=-(1.-QQ(1))*(1.+PQ(1IP))/Awt+(1./Bwt+1./Cwt)* 
$ PQCIP)*QQ( IP) 


CONTINUE 
CALL DIAGEQ(I2) 
FINOS CORRECTIONS 


*PQ(T)-(UCIP)-VCIP)+WwC(IP))*OQ(T) 


ADJU 
ADJU 
ADJYU 
ADJU 
ADJ 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADUJUU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADUJU 
ADUJU 
ADJU 
ADJU 
ADJU 
ADJU 


*- 8 © © 


60 


AW=1./SP(KOUNT, 1) 
BW=1./SDO(KOUNT, 1) 
CW=1./STCKOUNT, 1) 
UC(1)=SORT(ABS(U(1)+xX(1)*(1.+PQ(1))/AW)) 
VC(1)=SORT(ABS(VE1)-X(4)*PQ(1)/Bw)) 
WC(1)=SORT(ABS(W(1)4+X(1)*PO(1)/CW) ) 
AW=1./SP(KOUNT .N) 
BW=1./SD(KOUNT,N) 
CW=1./ST(KOUNT,N) 
UC(N)=SORT(ABS(U(N)-X(12)*( 1. -QQ(NM))/AW)) 
VC(N)=SQRT(ABS(V(N)-X(12)*QQ(NM)/BW) ) 
WC(N)=SORT(ABS(W(N)+X(12)*QQ(NM)/CW) ) 
I2=1 
DO 40 I=UJU,NM 
I2=12+1 
12M=I2-1 
AW=1./SP(KOUNT,1) 
BW=1./SD(KOUNT,I) 
CW=1./ST(KOUNT,I) 
IM=I-1 
IF(IT.E€Q.JJ)IM= 1 
UC(I)=ABS(UCT) +(-X(12M)*(1.-OOC IM) )+X(12)*(1.+PQ(1)))/AW) 
UC(1)=SQRT(UC(I)) 
VC(I)=ABS(V(T) - (XC 12M) *QQ(IM)+xX(12)*PQ(1))/BW) 
VC(I)=SQRT(VC(T)) 
WC(I)=ABS(W(I) +(X( 12M) *QQ(IM)+xX(12)*PQ(1))/cw) 
WC(I)=SORT(WC(I)) 
GETS ADJUSTED VALUES 
AQUUSTS ON TRIANGLE INEQUALITIES 
K=O 
DO 68 I=1,.N 
IF(I.GT.1.AND.I.LT.UJ) GO TO 68 
AU=UC(TI) 
AV=VC(1) 
AM=WC(I) 
AMAX=AMAX1(AU,AV,AM) 
EE=E1*AMAX 
EF=E2*AMAX 
AW=SP(KOUNT,I) 
BW=SD(KOUNT,1I) 
CW=ST(KOUNT, I) 
COR=AU+AV-AM-EE 
DIV=AW+BW+CW 
IF(COR.GT.O.) GO TO 60 
COR=(AU+AV-AM-EF)/DIV 
AU=AU~COR*AW 
AV=AV-COR*BW 
AM=AM+COR*CW 
GO TO 64 
COR=AU-AV+AM-EE 
IF(COR.GT.O.) GO TO 62 
COR=(AU-AV+AM-EF)/DIV 
AU=AU-COR*AW 
AV=AV+COR*BW 


ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
AODJU 
ADWJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJ 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
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ADJL 
ADWJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 


128 


62 


64 


66 


68 


100 


110 


112 
114 


1010 


AM=AM-COR*CW 

GO TO 64 

COR=-AU+AV+AM-EE 

IF(COR.GT.O.) GO TO 66 

COR=( -AU+AV+AM-EF)/DIV 
AU=AU+COR*AW 

AV=AV-COR*BW 

AM=AM-CORtTCW 

K=K+ 4 

UC(1)=AU 

VC(1I)=AV 

WC( I )=AM 

CONTINUE 

KMAX=K 
IF((ITER.EQ.0).OR.(KMAX.NE.O)) GO TO 110 
GO To 112 

ITER=ITER+1 

IF(ITER.LE.MAXIT) GO TO 12 

IF (ISS.NE.1) GO TO 999 

ITER=1 

ISS=2 

VTA=VC(1) 

WTA=WC(1) 

DO 120 I=JJ,NM 

IM=I-1 

IF(I.EQ.JJ)IM=1 

VIB=VC(I) 

WTB=WC(I) 

VC(L)=(VC(14+1)4+2. *VTB+VTA)*0.25 
WC(L)=(WC(14+1)+2. *WTB+WTA)*O0. 25 
VTA=VTB 

WTA=WTB 

CONTINUE 

GO TO 12 

CALCULATE THE CORRECTED VARIANCES 
DO 1010 I=1,N 
IF(I.GT.1.AND.1.LT.uuJ) GO TO 1010 
SP(KOUNT,I)=UC(I)**2 
SD(KOUNT, 1) =VC(1)**2 
ST(KOUNT,1I)=wc(1)**2 

CONTINUE 

RETURN 

END 

SUBROUTINE CHECK 
COMMON/CHK/P(4,4,3),RHO(4,4,3),NO(2) 
COMMON /WINCOM/OGH, FCORY ,DX5,0Y5,DUMMY( 14) 
COMMON /CHIC/LA(16),NB(2), IWSYM,USH, VSH, DUSH,DVSH 


NB(1) = O 
NB{(2) = 0 
CALL GROUP 
NS=0 
NR=1 


1F(NO(1).EQ.0.AND.NO(2).EQ.0) GO TO 1000 
DO 640 KL=1,2 


ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJ 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
AD JU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
ADJU 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
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130 
131 
132 
133 
134 
135 
136 
137 
138 
139 
140 
1414 
142 
143 
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
1614 
162 
163 
164 
165 
166 
167 
168 
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450 


SOO 


510 
520 


530 
540 


640 


1000 


+ © = 28 «© 


IF (NO(KL).EQ.0) GO TO 640 
CONTINUE 

NNR=4*NR 

IF(NO(KL).LE.NNR) GO TO 500 
NR=NR+1 

GO TO 450 

CONTINUE 

11=NR 

Jt=NO(KL)-CNR-1)*4 

SH1 = 6. 

SH2 = 6. 

DP = P(T1,J4,2) - P(T1,JU1,1) 
IF (DP) 510,520,510 

SH1t = ABS(P(T1,J1,2)/0P) 

DP = P{11,JU1t,2) - P(114,J1,3) 
IF (OP) 530,540,530 

SH2 = ABS(P(1I1,JU1,2)/DP) 
IF(SH1.LT.4.0.0R.SH2.LT.4.0) GO TO 640 
IF(SH1.GT.9.0.0R.SH2.GT.9.0) GO TO 640 
NR=4 

NS=NS+#¢ 

CONTINUE 

RETURN 

IWSYM = ICHAR(’*7) 

RETURN 

END 


SUBROUT INE CORLAT(A,B,C,D,E,F,G,H,AI,.AU,AK,SP1,SP2,SD1,S$D2,ST1, 
1 ST2,SU1,SU2,SV1,SV2,UD1,UD2,VD1,V02,RD,RT.RV) 
IF(SD1*ST1*SD2*ST2*RDY*RTtRV.GT.0.)GO TO 5 


DEFAULT VALUES AVOID DIVISION BY ZERO 
IF(SD1.LE.0.) SD1=0.001 

LECST 1. LEO. ST 102004 

IF(SD2.LE.0.) SD2=0.001 

IF(ST2.LE.0.) ST2=0.001 

IF(RD.LE.O ) RD = .0O001 


IF(RT.LE.O.) RT = .CQD001 
IF(RV.LE.O.) RV = .QO001 

CONTINUE 

IF (ABS(UD1).LE.0.) UD1 = 0.001 

IF (ABS(VD1).LE.0.) vO1 = 0.001 

IF (ABS(SU1).LE.O.) SU1 = 0.001 

IF (ABS(SV1).LE.0.) SV4 = 0.001 

IF (ABS(UD1).GE.1.) UD1 = 0.99*UD1/ABS(UD1) 
IF (ABS(VD1).GE.14.) VD1 = 0.99*VD1/ABS(VD1) 


A=RD*SD2/SD1 
B=SD2*SORT(1-RD*RD) 


TD2=(SP2*SP2-SD2*SD2-ST2*ST2)/(2*SD2*ST2) 
TDO1=(SP4*SP1-SN1*SD1-ST1*ST1)/(2*SD1*ST1) 


IF(ABS(TO1) LE.O.) TOt=.001 
IF(ABS(TO2) LE.O.) TO2=.001 


IF(ABS(TD2).GE.1.0) TO2=0.99*TD2/ABS(TD2) 
IFCABS(TD1).GE.1.0) TD1=0.99*TD1/ABS(TD1) 
C=(ST2/ST1)*(RT-RD*TD2*T01)/(1-TO1*TDt*RD*RD) 


D=(RT*ST2-CYST1)/(A*TD1I*SD1) 


CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CHEC 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORI. 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
CORL 
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10 


25 


10 
15 


20 


30 


E= ST2*ST2-C*C*ST1*5T1-D*D*SD2*SD2-2*C*D*RD*TDI*ST1*S02 CORL 
IF(E.GE.O.) GO TO 10 CORL 
E=0. CORL 
E=SORT(E) CORL 
F=(SU2/SU1)*(RV-RD*UD2*UD1)/(1-RO*RD*UD 1*U0 1) CORL 
G=(RV*SU2-F*SU1)/(RD*UD1*SD2) CORL 
H= SU2*SU2-F *F*SU1*SU1-G*G*SD2*SD2-2*F *G*RD*UD 1*SD2*SU1 CORL 
IF(H.GE.O.) GO TO 15 CORL 
H=O. CORL 
H=SQRT(H) CORL 
AI=(SV2/SV1)*(RV-RD*VD2*VD1)/(1-RD*RD*VD1*VD1) CORL 
AJ=(RV*SV2-AIlI*SV1)/(RD*VD1*SD2) CORL 
AK= SV2*SV2-AI*AL*SV1*SV1-AU*AU*SD2*SD2-2*AI *AU*RD*VD 1*SD2*SVICORL 
IF(AK.GE.O0.) GO TO 25 CORL 
AK=O. CORL 
AK=SQRT( AK) CORL 
RETURN CORL 
END CORL 
FUNCTION CORREL(X) CORR 
DATA A,B/19.51615854, 1.000416939/ CORR 
RHO = 1./EXP(B*tX) CORR 
IF(X .LT. 0.05) RHO = t. - A*X**2 CORR 
CORREL = RHO CORR 
RETURN CORR 
END CORR 
SUBROUTINE DIAGEQ(N) DIAG 
A(L,U)=DIAG. TERMS, I=ROW NO., J=DIAG. NO. DIAG 
B(I)=RIGHT SIDE TERMS DIAG 
N=NO. OF ROWS DIAG 
K=NO. OF BORDER DIAGONALS, M=K+1=INDEX OF PRIN. DIAG DIAG 
Q2KH=TOTAL NO. OF OIAGS. DIAG 
X(I)=SOLUTION DIAG 
COMMON/ADUCOM/A(26,3),B(26),X(26),KOUNT DIAG 
K= 1 DIAG 
M=K+ 4 DIAG 
DO 30 L=1,N DIAG 
ALM=A(L,M) DIAG 
A(L,M)=1. DIAG 
IF(L.EQ.N) GO TO 15 DIAG 
I12=MINO(K,N-L) DIAG 
pO 10 I=1,12 OIAG 
MP I =M+I DIAG 
A(L,MPI)=A(L,MPI)/ALM DIAG 
B(L)=B(L)/ALM DIAG 
IF(L.EQ.N) GO TO 30 DIAG 
Do 25 I=1,12 DIAG 
LPI=L4I DIAG 
FACT=A(LPI,M-I1) DIAG 
DO 20 J=1,12 DIAG 
MJI=M+tJ-I DIAG 
A(LPI,MJI)=A(LPI,MJI )-ACL, M+) *FACT DIAG 
B(LPI)=B(LPI)-B(L)*FACT DIAG 
CONTINUE DIAG 
X(N) =B(N) DIAG 


10 


40 
90 


NM1=N- 1 
DO 50 L=1,NM1 
NML =N-L 
SUM=O. 
I2=MINO(K,L) 
DO 40 I=1,12 
SUM=SUM+A(NML ,M+I)*X(NML+I) 
X (NML )=B(NML )-SUM 
RETURN 
END 
SUBROUTINE FAIR(PG,0G,TG,PU,DJU,TJ,IH,P,D.T,DPXG, 


& DPYG.DPXJU,DPYJU,DPX, DPY .DTXG,DTYG,DTXJU,DTYJU,DTX ,DTY) 


FAIRS BETWEEN ZONAL MEANS AND JACCHIA VALUES 90 LE HEIGHT LE 
120 KM 
DIMENSION CZ(7) 
FAIRING VALUES 
DATA CZ /1.0,0.9330127,0.75,0.5,0.25,0.0669873 ,.0.0/ 
HEIGHT INDEX 

I =(IH - 85)/5 
ZONAL MEAN FAIRING COEFFICIENT 
CZI = C21) 
JACCHIA FAIRING COEFFICIENT 
SZI1 = 1.0 - CZI 
FAIRED TEMPERATURE 
T = TG*cZI + TUtS§Z!I 
FAIRED DENSITY 
D = EXP(ALOG(DG)*CZI + ALOG(DJU)*S7Z1) 
FAIRED GAS CONSTANT 
RG=PG/(0G*TG) 
RU=PU/(DuUtTY) 
R=( CZI*RG ) + ( SZI * RJ ) 
P= R*D+tT 

DPX = DPXG*CZI + DPXUJ*SZI 
DP/DY FOR GEOSTROPHIC WINDS 
OPY=DPYG*CZI+DPYuU*SZI 

DTX = DIXGtCZI + DTXJU*SZI 

DT/DY FOR THERMAL WINDS 
DTY = OTYG * CZ1 + DOTYJ * §Z!I 
RETURN 
END 


SUBROUTINE FAIRS ( DHEL1 ,DHEL2 ,DLG1 ,DLG2 ,IH ,FDHEL ,FDLG ) 


DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
DIAG 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FATR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAIR 
FAER 
FAIR 
FAIR 
FARS 
FARS 
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C** THIS SUBROUTINE FAIRS BETWEEN THE REGION ABOVE 500 KM, WHICH 


C** INVOKES THE SEASONAL-LATITUDINAL VARIATION OF THE HELIUM NUMBER 


C** DENSITY (SUBROUTINE SLVH ), AND THE REGION BELOW, WHICH DOES 
C** NOT INVOKE ANY SEASONAL-LATITUDINAL VARIATION AT ALL. 


Cr 
C** INPUTS: DHELt = HELIUM NUMBER DENSITY BEFORE INVOKING SLVH 
Chee DHEL2 = HELIUM NUMBER DENSITY AFTER INVOKING SLVH 
cnt DLG1 = TOTAL DENSITY BEFORE INVOKING SLVH 

Ce DLG2 = TOTAL DENSITY AFTER INVOKING SLVH 

ost TH = HEIGHT ( KM ) -- INTEGER 

Gt IBFH = BASE FAIRING HEIGHT ( KM ) -- INTEGER 


**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 
**FARS 


ctr 
C** 
C** 


OUTPUTS: 


FOHEL 
FOLG 


= FAIRED HELIUM NUMBER DENSITY 
= FAIRED TOTAL DENSITY 


++ FARS 
**FARS 
**FARS 


Cee RE OR RRR RK ET ORR ERY KF RR EKER ECR ER RRR EH KEE KEPARS 


PARAMETER (IBFH = 
DIMENSION CZ ( 6 ) 
DATA CZ / 1.0, ©0.9045085, 0.6545085, 0.3454915, 0.0954915, 0.0 /FARS 


HEIGHT INDEX 


I = { 


CZI = 


IH - 


om came | 


440) 


IBFH ) /10 + 1 
NON-SLVH FAIRING COEFFICIFNT 


I) 


SLVH FAIRING COEFFICIENT 


SZI = 


120 


FAIRED DENSITY 


FOLG 


= { DLGt 


CZI 


FATRED HELIUM NUMBER DENSITY 
) + ( DHEL? * SZI 


FOHEL 


RETUR 
END 


= ( DHEL1 


N 


* CZi 


‘ CZI ) + ( DLG2 


¥ 


S71 ) 


SUBROUTINE GAUSS(Z2,NMIN,R,TX,.T1,73,T4,A2) 


) 


FARS 
FARS 
FARS 


FARS 
FARS 
FARS 
FARS 
FARS 
FARS 
FARS 
FARS5 
FARS 
FARS 
FARS 
FARS 
FARS 
FARS 
GAUS 
GAUS 


COI II IN IORI IORI HORROR IIIT IORI I TOI TOT TE TIE EE EEE ET ETT RARE T HERE EGAUS 


SUBDIVIDE TOTAL INTEGRATION-ALTITUDE RANGE INTO INTERVALS 
SUITABLE FOR APPLYING GAUSSIAN QUADRATURE , 


Cc** 
Cr 


C+ 
C*+ 


C 


C 


OF POINTS FOR INTEGRATIONRFOR EACH SUB- INTERVAL 


PERFORM 


GAUSSIAN QUADRATURE. 


CR OO HO OE TR RT ET THK EH HHT TE EGAUS 


. C18,6). 
NGAUSS , 


X(8,6), 
NMIN 


<u 


ALTMIN/S90.,105.,125.,160.,200. .300 


REAL ALTMIN (9) 
INTEGER NG (8) , 
DATA 

DATA NG / 4 


ee, 


. 6 


’ 


COEFFICIENTS FOR GAUSSIAN QUADRATURE 


DATA 


Cc / 


,5955906 


. COO0000 
. 314178548 
. OO00000 
. 2369269 
. 2369269 
.1713245 
. 3607616 
. 1294850 
.3818301 
. 1012285 
. 3626838 


.8888889 
.OO00000 
.6521452 
. 0000000 
.4786287 
. OOG0000 
. 3607616 
.1713245 
. 2797054 
.2797054 
. 2223810 
. 3137067 


ABSCISSAS FOR GAUSSIAN QUADRATURE 


6 , 6 


(9555556 
.OO00000 
.652 1452 
. OO00000 
.5688889 
. OOO0000 
.4679139 
. OOO0000 
. 3818301 
. 1294850 
.3137067 
. 2223810 


MOLWE 


SET THE NUMBER 


~ AND THEN 


., 900. , 1500 


6 ., 6 


. 0000000 
. 0000000 
.3478548 
,OOO0000 
.4786287 
. OCOOO000 
.4679139 
. OOOO000 
.44179592 
. COO0000 
. 3626838 
. 41012285 


_,2500./ 
6 / 


**GAUS 
**GAUS 
**GAUS 
**GAUS 


GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 


11 


12 


DATA X / -. 
. CO00000 
.8611363 
. CO00000 
.906 17938 
.906 1798 
.9324635 
.6612094 
.9491079 
.4058452 
.9602899 
. 1834346 


7745967 


. OOO0000 
. OO00000 
.3399810 
. 0000000 
.5384693 
. 0000000 
.6612094 
.9324695 
. 7415312 
.7415312 
. 7966665 
.5255924 


.7745967 
. 9000000 
. 3399810 
. 0000000 
. OOO00000 
0000000 
. 2386192 
. OOO0000 
.4058452 
.9491079 
.5255324 
. 7966665 


. OO00000 
. 0000000 
.8611963 
- OO00000 
.5384693 
. 0000000 
. 2386192 
. OOO0000 
. QOO0000 
. DO00000 
. 1834346 
. 9602899 


GRAVITY(ALT)=9.80665/((1.+ALT/6.356766E3 )**2 ) 


R = 0.0 
DO 2 K = NMIN , 8 
NGAUSS = NG (K) 
A = ALTMIN (K) 
D = AMIN1 ( Z2 , ALTMIN (K+1) ) 
RR = 0.0 
DEL = 0.5 * (D0 - A ) 
J = NGAUSS - 2 


DO i I = 1 , NGAUSS 


Z = DEL * ( X€I,U) + 1. ) + A 
RR=RR+C(1,JU)*MOLWT(Z)*GRAVITY(Z)/TEMP(Z,TX,T1,73,T4,.A2) 


CONT INUE 
RR = DEL * RR 
R = R + RR 
IF ( D .EQ. 22 ) RETURN 
CONTINUE 
RETURN 
END 


SUBROUTINE GEN4D 
Maree GENERATES NG = 9 OR 16 4D PROFILES P,0,T AND SIGMAS SP,SO,ST AT 
GRID OF LATITUDES ANDO LONGITUDES GLAT,GLON. CURRENT LATITUDE, 
LONGITUDE=CLAT,CLON. PREVIOUS LATITUDE,LONGITUDE=PLAT,PLON. 
COMMON/C4/GLAT(16),GLON(16),NG,P(16,26) ,.0(16,26),T( 16,26), 
$ SP(16,26),S5D(16,26),ST(16,26),PLON,CLON,HS 
COMMON/ IOTEMP/TOTEM1, IOTEM2,1IUS,DDD, XMJD,PLAT,CLAT, 
$ NSAME ,RP1,RD1.RT1,5P1,.501,ST1,RU1,RV1,SU1,SV1, 
$ MN, IDA, ITYR,H1,PHIIR, THETIR,GZ,RI,Z2,PHIR, THETR,F10O,F10B8,AP, 
$ IHR,MIN,NMORE ,DX,HL,VL,.02,8,EPS,IOPP,LOOK,DUMMY(20) 
COMMON /PDTCOM/1U4,MONTH,IOPR, 
PSP(15,19,18),0SP(15, 19,18),TSP(15,19,18),USP( 15,19, 18), 
vSP(15,19,18), 
PG(21,19),0G6(21,19),TG(21,.19),UG(21,19). 


GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GAUS 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
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+ 
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1OO 
110 


120 


130 


PAQLI7,.5), DAQ( 17.5), 1TAQ9017.5),UAQ(17.5),VAQ(17.5), 
poo(t7.5).poec 17,5), 1T0eC17,5),U00( 17,5), VDO(17,5). 
PP(29. 19), PP(29,19),1R(29, 19), UR(29, 19), VR(29,19), 
Py, DO,1Q,UO,VG,POA, DOA, TOA,UA,VA,1OPQ 
PLP(25, 10), DLP (25, 10), 1LP025, 10), ULP(25, 10), VLP(25,10),UDL(25, 
19), VOL (25, 19).UDS(25. 10), VbS€25, 10) 
COMMON! /ADUCOM ‘DUM( 190) .KOUNT 
COMMON/TFRIP, IPRE,NLIMET 
DIMENSION NGONN( 26) 
IF (NSAME FQ 1) RETURN 


LOOK-© 

Fo O.O87453929923 
NG = 16 

Dx = FION - CION 


IF (OX. GT. 190. DX = OX - 3GO. 

IF COX LE. t8O. OX = DX 4 360. 
LONGITUDE DISPLACEMENS FROM PREVIOUS TO CURRENT POSITION 
DY = CLAT - PLAT 
LATITUDE DISPLACEMENT FROM PREVIOUS TO CURRENT POSITION 
TE CBi4 2 20.10.20 
jf eee ee Ws Os are am 


K = O 

GD TO 10 

THETA > 1890. + STIGN(90. .DxX) 
GON TQ 30 

THETA = ATAN(DX/DY)/F 


IF (DY.GT.O0.) JHETA = THETA + 180. 

IF €(THETA.LT.O.) THETA = THETA + 360. 

THETA = AZIMUTH ANGLE OF TRAJECTORY, USED TO ORIENT LAT-LON GRID 
K = INV((THETA # 67.5)/45.) 

INDEX USFD If! COMPLITED GO TO FOR 110 THRU 180 

IF (K.G1.8) KeK-B 

NORTH Fol AR GRID 

IF (CLAT.G1.75.Q0}G0 TO 200 

SOUTH POLAR GRID 

IF (CLAT.LI.-75.9)G9 Ta 200 


INITIAL ESTIMATE OF REFERENCE LATITUDE (LOWER LEFT GRID POINT ) 


LATO - SY*INT(CLAT/S. ) 
IF (CLAT.LE.O.) LATO = LATO - 5 


-TINETIAL FSPTIMATE OF RFEFFRENCE LONGITUDE (LOWER LEFT GRID POINT) 


LONO-S'IHT cet ons. ) 

ADWIUSTS LATO.LONO ACCOPDING TO DIRECTION OF TRAJECTORY AZIMUTH 
If (K GT.O) AN TO 100 

LATO - LATO 5 

LONO- LONO + 10 

co TO 190 

GO TO £110, 120, 130, 149, 150,160,170, 180),K 
LALO. SO LAIOMO 

LONOC - LONG + 10 

GQ 10 t990 

LATO = LATO-10 

LONO - LONGO? 15 

GO 10 190 

LATO - LATO-5 


GEN4 
GENA 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN-+4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GENA 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GENA 
GEN4 
GEN4A 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 
GEN4 


13 


14 


140 


150 


160 


170 


180 


190 


192 


. 5 6 @ 8 


.POLAR GRID LATITUDES 1-8 = +75 (N)} OR 


. = 8 & 4 


325 
330 


LONO = LONO+15 
GO TO 190 

LONO = LONO+15 

GO 10 190 

LONO = LONO+10 

GO TO 190 

LONO = LONO+5 

GO TO 190 

LATO = LATO-5 

LONO = LONO+5 

GO To 190 

LATO = LATO-10 

LONO = LONO+S5 
*TF(LONO.GE.360)LONO = LONO - 360 
IF(LATO.GT.75) LATO=75 

DLI=1.25 

IF(ABS(CLAT).GE.18) GO TO 192 
DLI=3.0 


LATO=- 18 
DO 195 I=1,4 
I12 = [+12 


DO 195 J=1I,112,.4 
GLAT(JU)=LATO+DLI *(JU-I1) 


LATITUDE, LONGITUDE GRID AT 5 DEGREE INTERVALS 


GLON(J) = LONO - 5. * (I - 1) 
GO TO 400 

POLAR GRID 

NG = 9 

DO 210 J=1.8 


GLAT(JU) = SIGN(75.,CLAT) 


POLAR GRID LONGITUDES 1-8 AT 45 DEG INTERVALS 


GLON(J) = 45.*(uJ-1) 


POLAR GRID LATITUDE 9 = POLE +90 OR 


GLAT(9) = SIGN(90. ,CLAT) 
POLAR GRID LONGITUDE 9 = O 
GLON(9) = O. 


GENERATES 16 PROFILES (OR 9 PROFILES FOR POLAR GRID) 


CALL GRID4D 
DO 390 I = 1 
DO 330 JU = 1 
NGOOD(J) = 1 


.NG 
26 


IF(P(I,J).LE.0.0.0R.D(1I,JU).LE.0.0.0OR T(I,JU).LE.0.0)NGOOD(uU)=0 


IF(NGOOD(J).EQ.0)GOTO 330 
RATIO = PCI,JU)/(D(I.U)*T(I,d)) 


IF(RATIO.GT.286.0.AND.RATIO.LT.288.0}GOTO 330 


NGOOD(J) = O 


IFCIPRT.EQ.O)WRITE(6,325)I,.U,RATIO 
FORMAT(’ GAS LAW VIOLATION. I.U,RATIO 


CONTINUE 
DO 340 J = 3,26 


IF (NGOOD(JU).EQ.0.0R.NGOOD(J-1).EQ.0)GOTO 340 


DENOM = 1./T(I,WU) 


IF(CABS(T(CI,J)-TC1,J-1)).GT.0.01 )DENOM=ALOG(T(1I,JU-1)/T(I.J))/ 


‘,214,G12.4) 
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& (T(I,J-1)-T(I,J)) 
RATIO = ALOG(P(I,JU-1)/P(1I,U))/DENOM 
IF(RATIO.GT.30.7.AND.RATIO.LT.37.6)GOTO 340 
NGOOD(JU) = O 
IFCIPRT.EQ.O)WRITE(6,335)1,JU,RATIO 
335 FORMAT(’ HYDROSTATIC VIOLATION. I,JU,RATIO = ',214,G612.4) 
IF(JU.EQ.26)GO TO 345 
K1=J5+ 1 
DO 336 K 
336 NGOOD(K) 
GO TO 345 
340 CONT INUE 
345 NBAD = O 
DO 360 JU = 1,26 
IF(NGOOD(J).GT.O)GOTO 360 
NBAD = NBAD + 1 


K1,26 
0 


P(I,J) = O. 
D(I,J) = O. 
T(I,JU) = O. 


360 CONTINUE 
IF(NBAD.LE.NLIMIT) GOTO 390 
WRITE(6,380) 
38O FORMAT(’ UNABLE TO GENERATE 4-0 GRID. TOO MANY TEST VIOLATIONS’ ) 
STOP 
390 CONTINUE 
DO 600 I=1,NG 
IHV O 
SPR 0.0004 
SDR 0.0004 
STR 0.0004 
DO 420 JU = 8,26 
CHECK = 1. 
IF(P(1,JU).LE.0.0.0R.SP(1I,JU).LE.0.0)CHECK 
IF(D(I,J).LE.0.0.OR.SO0(I,JU).LE.O0.O0)CHECK O. 
IF(T(I,J).LE.O.0.OR.ST(I,U).LE.0.0)CHECK O. 
Cee sae FINDS INDEX IHV OF LAST HEIGHT ABOVE 6 KM WITH NON-ZERO DATA 
IF(CHECK.GT.O.)GO TO 420 
IHV = J-1 
GO TO 440 
420 CONTINUE 
C HEIGHT = HEIGHT INDEX - 1 
440 21 = IHV -1, 
IFC IHV.EQ.0)GO TO 491° 
Cc SPR,SDR,STR=SIGMAS AT HEIGHT 21 
SPR = SP(I,IHV} 
SDR=SD(I, IHV) 
STR=ST(I,IHV) 
IF(SPR.LE.O.O)SPR 
IF(SDR.LE.O.0)SOR 0.0004 
IF(STR.LE.O.O)STR 0.0004 
IF(IHV.GT.26-NLIMIT)GOTO 441 
WRITE(6,442)IHV 
442 FORMAT(’ UNABLE TO GENERATE 4-D GRID. IHV = ’,19) 
STOP 
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Ad{ CONTINUE GEN4 
Sah ane IF HEIGHT Z1 GEQ 20 KM, USE ZONAL MEAN AT 30 KM FOR INTERPOLATION, GEN4 
C OTHERWISE USE ZONAL MEAN AT 25 KM GEN4 

IF (IHV.GE.21) GO TO 480 GEN4 

C.... EVALUATES ZONAL MEANS AT 25 KM FOR INTERPOLATION AND GEN4 
Cc FILL IN OF ZERO DATA GEN4 
CALL GTERP(25,GLAT(1I).P2,02,T2,PG,D0G,TG.DPY ,OTY,U2,UG) GEN4 

IHP = IHV + 1 GEN4 

DO 450 K=IHP,26 GEN4 

A AVOIDS INTERPOLATION OF P,D,T IF ONLY SIGMAS ARE ZERO GEN4 
IF(P(I,K).GT.0.0.AND.0(1,K).GT.0.0.AND.T(1I,K).GT.0.0)GO TO 445 GEN4 

H=K-4 | GEN4 

C.... INTERPOLATES BETWEEN 4D AT HEIGHT Z1 AND ZONAL MEAN AT 25 TO FILL GEN4 
C IN MISSING DATA GEN4 
CALL INTER2(P(1,IHV),DCI,IHV), TCI, IHV),2Z1,P2,02,7T2,25..PH,DH, TH,H)GEN4 
IF(IPRT.EQ.O)WRITE(6,583)I1,«K,PH,DH, TH GEN4 
P(I,K)=PH GEN4 
D(1.«)=DH GEN4 
T(I,K)=TH GEN4 

445 SP(1,K) = SPR GEN4 
SD(I,K)=SDR GEN4 
C.....SETS MISSING SIGMAS EQUAL TO SIGMAS AT HEIGHT 21 GEN4 
450 ST(I,K)=STR GEN4 
GO TO 491 GEN4 

C.... EVALUATES ZONAL MEANS AT 30 KM FOR INTERPRETATION AND FILL IN OF’ GEN4 
C ZERO DATA GEN4 

480 CALL GTERP(30,GLAT(1I),P2,02,.T2,PG,0G,TG,DPY ,DTY,U2,UG) GEN4 

C COMPUTE PERTURBATIONS TO ZONAL MEAN MODEL GEN4 

CALL PDTUV(PSP,DSP,TSP,USP,VSP,.GLAT(I),GLON(I).30.0P,D0O,DT, GEN4 

$ DPX ,DPY ,DTX,DTY,DU.DV) GEN4 
C.... ADD STATIONARY PERTURBATIONS TO ZONAL MEAN MODEL GEN4 
P2 = P2*(1. + DP) GEN4 

D2 = D2*(1. + OD) GEN4 

T2 = T2*(1. + OT) GEN4 

IHP = IHV + 1 GEN4 

DO 490 K=IHP,26 GEN4 
C.....AVOIDS INTERPOLATING P,D,T IF ONLY SIGMAS ARE ZERO GEN4 
IF(P(1,K).GT.0.0.AND.D(1I,K).GT.0.0.AND.T(I,K).GT.0.0)GO TO 485 GEN4 

H=K-1 GEN4 
Ces INTERPOLATES BETWEEN 4D AT HEIGHT 71 AND GROVES AT 30 KM TO GEN4 
C FILL IN MISSING DATA GEN4 
CALL INTER2(PC(I,IHV),O(I,IHV),. TCI. IHV),21,P2,02,72,30. ,PH,DH, TH.H)GEN4 

IFC IPRT.EQ.O)WRITE(6,583)1,K,PH,OH,TH GEN4 
P(I,K)=PH GEN4 
D(I,K)=0H GEN4 
T(I1,K)=TH GEN4 

485 SP(I,K) = SPR GEN4 
SD(I,K)=SDR GEN4 

C SET MISSING SIGMAS AT HEIGHT 1 GEN4 
490 ST(1,K) = STR GEN4 
491 CONTINUE GEN4 
IHP = IHV - 1 GEN4 

SPO = SP(I,1) GEN4 

S00 = SDfI,1) GEN4 
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STO = ST(I,1) 


IF(SPO.LE.O0.0)SPO = 0.0001 
IF(SDO.LE.0.0)SDO = 0.0001 
IF(STO.LE.O.0)STO = 0.0001 


DO 492 K = 1,9 
IF(SP(I,K) .LE. O.) SPC(I,K) SPO 
IF(SD(I,K) .LE. O.) SD(L,K) ~ SDO 
IF(STCI,K) .LE. O.) ST(I,K) = STO 

DO 495 K=10,I1HP 

SETS ALL ZERO SIGMAS TO SIGMA AT HEIGHT 21 


It 


IF (SP(I,K).LE.O.0.AND.P(I,K).GT.O.) SP(I,.K) = SPR 
IF (SO(I,K}.LE.0.0.AND.D(1,K).GT.O.) SO(I1,K) = SDR 
IF (STC(I,K).LE.O.O0.AND.T(I.K).GT.O.) STfI,K) = STR 


PA = P(I,1) 
TA = T(I,1) 


R =287.05 

G = GZ*(1.+(Z/(RI-Z)))**2 
K = 2 

PB = P(1,K) 


TB = T(I,K) 
IF ((PB*TB) .GT. 0.) GO TO 520 


K =KK + 1 
GO TO 510 
IF(ABS(TA-TB).LE.0.01)GOTO 570 

IF(TA*TB.LE.O.0)GO TO 570 

TZ = (TA-TB) / ALOG(TA/TB) 
GO TO 575 
TZ = TA 

HS = K - 1. 

IF(PB*PA.LE.0.0)GO TO 576 

HS = K - 1. + ©0.001*R*TZ*ALOG(PB/PA)/G 
KM = K - 2 


IF(ABS(K-1-HS).Gl.0.1) GO TO 578 
GAM=TB-T(I,K+1) 
IF(ABS(GAM).LE.0.01)GOTO 590 
GO TO 582 
IF(ABS(TA-TB).LE.0.01)GOTO 590 
GAM=(TA-1TB)/(K-1-HS) 
KM1=KM+ 1 
IF(ABS(GAM).G1.G) GAM=SIGN(G, GAM) 
DO 585 JD=1,KM1,1 
J=JD- 1 
TJU=TA-GAM* (J-HS ) 
PJU=PA*(TJ/TA)**(G/(R*GAMtO.001)) 
DU=PU/(R*TU) 
IFC IPRT.EQ.O)WRITE(6,583)1,JU+1,PU,Dd, TJ 
FORMAT(’ NEW VALUES COMPUTED AT I.J P.D,T ’,214,F9.0,F9.5,F9.2) 
P(I,J+1)=Pud 
D(I,JU+1)=Du 
TCI, Ju+14)=TJ 
GO TO 599 
KM 1=KM+ 1 
DO 595 JD=1,KM1, 1 
J=JD- 1 
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TU=TA 
PU=PA*EXP(-G*(JU-HS)/(R*O0.001*TU)) 
DJU=PU/(R*TUJ) 

IFC IPRT.EQ.O)WRITE(6,5839)1,U+1,Pu.0U,TU 
P(I,J+1)=PuJ 

O(1I,.J+1)=Dd 

TCT, J+ )=Td 

IF (NSAME .EQ.2) NSAME=1 

HS=O. 

KOUNT = I 

CALL ADUJUST 

CONT INUE 

RETURN 

END 

SUBROUTINE GETNMC 


READS OATA FILE FOR NMC GRID NUMBERS, 
AND WRITES SCRATCH FILE FOR USE BY SELEC4. 


DIMENSION IP(15) 
CHARACTER N*2, IDUMMY*2 


COMMON/IOTEMP/ IOTEM1, IOTEM2,IUS, JDUM(6O) 
NREC=0O 


READ(IUS,300, END=90) N,IP 

FORMAT(A2, 1517) 

IF(N.NE.‘N’) GO TO G6 

DO 4 I=1,15,3 

M=IP(I) 

IF(M.LT.1) GO TO 5 

TU=IP(1+1)* 1000+IP(1+2) 

WRITE(IOTEM2) TJ 

NREC=NREC+ 1 

CONT INUE 

GO To 1 

IF(NREC.NE.1977) GO TO 6 

RETURN 

WRITE(6,200) NREC, IOTEM2 

FORMAT(1H1/1X,.16.’ RECORDS WRITTEN BY GETNMC IN SCRATCH FILE’ ,13) 
STOP 

WRITE(6,400) TUS 

FORMAT(’1 PREMATURE END-OF-FILE FOUND ON UNIT ’,12/ 
1‘0 CALLED FROM SUBROUTINE GETNMC.’ ) 

STOP ’ GETNMC ERROR.’ 

END 

SUBROUTINE GRID4D 

REAL LAT,LON 
COMMON/C4/LAT(16),LON(16),NP,P(16,26).R(16,26),T(16,26).5P(16,26), 
& SR(16,26),ST(16,26),DU1,DU2,DUMMY 

COMMON /POTCOM/IT,MONTH, IDUM, DUMMY 1( 30851) ,DUM2(2261) 
COMMON /IPRTP/IPRI ,NLIMIT 
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SUBROUTINE TO SELECT PRESSURE, TEMPERATURE, AND DENSITY PROFILES (GRID 


TOGETHER WITH THE NORMALIZED VARIABLES IN FACHI 
AT LAT/LONS SELECTED BY CALLING PROGRAM. 


USES DISK FILES MADE FROM NASA HUNTSVILLE MSFC 4°0 DATA TAPES 


COMMON/ IOTEMP/IOTEM1, IOTEM2,OUMMY2(61) 
COMMON /POINT/ IPT(16,5),.LL(16),0xXY¥(16,2) 
COMMON /ORDER/ IPTN( 16,5), IREAD(65,3) 


COMMON /INT/ D(208,5),1G(5),D0YX(2),DLA(4),0L0(4) 


NAME OF 4-0 FILE TO OPEN 
CHARACTER*7 NAME4D, FILES( 12) 
REAL DTEMP(213) 


INITIALIZE 


DATA FILES/’DAT4D01’ , ’DAT4D02’ , ‘DAT4003’ , ‘DAT4D04' , ’DAT4D05’ , 
&’DAT4D06/ , ‘DAT4D07' , ‘DAT4008’ , ’DATADO9’ , ‘DAT4D10’, ‘DAT4DI1’, 
&‘DAT4D12’/ 


ZERO=0.0 
ONE=1.0 
TEN=10.0 
HUNDR=100.0 
THOU= 1000.0 


OPEN 4-D FILE FROM MONTH VALUE 


NAME 4D=FILES( MONTH) 


OPEN( IT, FILE=NAME4D, STATUS=’OLD’ ,ACCESS=’DIRECT’, 


1FORM= ‘FORMATTED’ ,RECL=1269) 


APPROPRIATE 4-D INPUT FILE NOW POSITIONED - STORE NEEDED PROFILES 


20 CALL SELEC4 


29 


19 
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DO 29 IRN=1,64 

IFC IREAD(IRN,3).EQ.0)GOTO 29 
II = IREAD(IRN, 1) 

JJ = IREADCIRN,2) 


IF(IPRT.EQ.O)WRITE(6, 169) IRN, (IREAO(IRN,J),J=1,3), 
& IPT(II,uJ), IPTN(II, Jd) 


CONTINUE 

DO 38 II 1,NP 

po 19 I = 1,208 

p(1,5) = 0.0 

DO 31 J =1,4 ~~. 

IF (IPTN(II,J) .EQ. O)GOTO 31 

IRN = IPTN(II,J) - 9000 
IFCIPRT.EQ.O)WRITE(6,169)II,JU,IRN 
FORMAT(’ GRID4D’ ,6I6) 


READ IN RECORDS OF POINTS NEAR DESIRED LAT/LON 
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GRID 64 


21 READ(IT,22,REC=IRN) DOTEMP GRID 65 
22 FORMAT(208F6.2,2F4.2,3F4.0) GRID 66 
IF((OTEMP(213).NE.MONTH).OR. (OTEMP(212) NE.IPT(II,u)})) GOTO 39 GRID 67 

DO 30 1=1,208 GRID 68 

30 D(I,J)=OTEMP(1) GRID 69 
DLA(J)=DTEMP(209)*TEN GRID 70 
DLO(JU}=DTEMP(210)*TEN GRID 71 

31 CONTINUE GRID 72 
GRID 73 

IFCIPRT .EQ.O)WRITE(G6,99)II,LATC(II),LON(II),LLCII), CIPTCIE,Y), GRID 74 

& DLA(J),OLO(JU),0(0 1,9) *HUNDR ,D( 157,JU)/THOU,D(53,J),JU=1,4) GRID 75 
99 FORMAT(’ 4-D DATA ’,13,2F7.2,110/4(15,2F7.2,F8.0,F8.5,F8.2/)) GRID 76 
IF NECESSARY, INTERPOLATE GRID 77 
GRID 78 

LALO=LL(11) GRID 79 

DO 33 1=1,5 GRID 80 
IG(I)=IPT(II,1) GRID 81 

33 CONTINUE GRID 82 
IF(IG(2).NE.0) GO TO 35 GRID 83 

DO 34 I=1,208 GRID 84 
D(1,5)-D(1,1) GRID 85 

34 CONTINUE GRID 86 
GO TO 37 GRID 87 

35 IF(IG(5).NE.2) GO TO 36 GRID 88 
DYX(1)=DOXY(II,1) GRID 89 
DYX(2)=OxXY(II,2) 3 GRID 30 
GRID 91 

36 CALL INTRP4 (LALO) GRID 92 
GRID 93 

37 DO 38 I=1,26 GRID 394 
P(1I,1)=D(1,5)*HUNDR GRID 95 
R(II,1)=D(1+156,5)/THOU GRID 96 
TCELEI,E) =O0€1452,5) GRID 97 
DIVIDE=ONE GRID 98 
IF(P(II,1).GT.ZERO) DIVIDE=(PC(II,1)/HUNDR)**2 GRID 99 
SP(ITI.1)=D0(1+26,5)/DIVIDE GRID 100 
SP(IE,IT)=ABS(SP(II,1)) GRID 101 
DIVIDE=ONE GRID 102 
IF(R(II,I).GT.ZERO) DIVIDE=(THOUTR(II,1))**2 GRID 103 
SR(II,1})20(1+182,5)/DIVIDE GRID 104 
SRC II,1)=ABS(SR(II,1)) GRIO 105 
DIVIDE =ONE GRID 106 
IF(T(II,1).GT.ZERO) DIVIDE=T(IT,1)**2 GRID 107 
ST(II,1)=O(014+78,5)/DIVIDE GRID 108 
ST(CII,I)=ABS(ST(II,1)) GRID 109 

38 CONTINUE GRID 110 
IF(IPRT .NE.O)RETURN GRID 111 
WRITE(6, 2000) (LAT(1),.1=1,NP) GRID 112 
WRITE(6,2001)(LON(I),I=1,NP) GRID 1413 
WRITE(6, 2007) GRIC 144 

DO 501 I = 1,26 GRID 1415 

IH = I - 1 GRID 116 


WRITE(6, 2004 )IH, (SORT(SP(U,1}).JU=1.NP) GRIO 117 
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WRITE(6, 2002)IH, (P(JU,1),J=1,NP) 
WRITE(6, 2003) 

DO 502 I = 1,26 

IH =1- 1 

WRITE(6, 2004)IH, (SQRT(SR(U,1)),JU=1,NP) 
WRITE(6, 2004) IH, (R(J,1),JU=1,NP) 
WRITE(6, 2005) 

DO 503 I = 1,26 

IH =I- 1 

WRITE(6, 2004 )IH. (SQRT(ST(U,1)),U=1.NP) 
WRITE(6. 2006)IH, (T(U,I),JU=1,NP) 
FORMAT(’ 4-D DATA READ IN FROM FILE‘/’ LATITUDE’/3X, 16F8.3) 
FORMAT(’ LONGITUDE ’/3X, 16F8.3) 
FORMAT(’ PRESSURE’ ) 
FORMAT(1X,12,16F8.0) 
FORMAT(’ DENSITY’ } 

FORMAT(1X,12, 16F8.5) 
FORMAT(’ TEMPERATURE’ ) 

FORMAT(1X,12, 16F8.2) 


RETURN 
WRITE(G,40) IT, IRN,DTEMP(213),MONTH,OTEMP(212),11T,J,IPT(II,Y> 
FORMAT(’ ***** UNIT NO. 7,13,’ 1N ERROR.’, 
’ IRN =’,15,’ OTEMP(213) =’,13,’ MONTH -’,13,’ OTEMP(212) 
f TPTC*.12,°. 7.11.) e715) 
STOP 
END 


SUBROUTINE GROUP 

DIMENSION KOU(2) 
COMMON/CHIC/LA(4,4),NB(2), IWSYM,USH, VSH,DUSH,DVSH 
COMMON /CHK/P(4,4,3),DEN(4,4,3).NO(2) 
COMMON/WINCOM/DGH,FCORY ,Ox5,DY5,0UMMY( 14) 
FCORX = FCORY*DX5/DY5 

KK=1 

DO 100 I1=1,4 

DO 100 J=1,4 

LACSI,J)=4*(1I-1t)+uJ 

CONTINUE 

CONTINUE 

DO 250 M=1,4 

DO 250 N=1,4 

IF (KK.EQ.1) GO TO 210 

I=5-M 

J=5-N 

NN=- 4 

N4=- 4 

GO TO 220 

CONT INUE 

I=M 

J=N 

NN= 1 

N4=1 

CONTINUE 

IF (N.EQ.4)}) GO TO 225 

DINX=FCORX* (DEN(CI,J+NN,2)+#+DEN(I,JU,2))/2 
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IF(ABS(DINX).LE.O.)GO TO 225 
VY=(P(I,J+NN, 2)-P(1,JU.2))/DINX 
IF(ABS(VY).GT.100) GO TO 225 
LA(1,J)=MINO(LA(I,J),LACI,J+NN) ) 
LAC I, J#NN)=LA(I,J) 

CONTINUE 
IF(M.EQ.4) GO TO 250 
DINY=FCORY* (DEN(CI4N4,JU,2)+DEN(I.JU,2))}/2 

IF(ABS(DINY).LE.O.)GO TO 250 
VX=(PCI4N4,U,2)-P(I,U,2))/DINY 
IF(ABS(VX).GT.100) GO TO 250 
LACT, J)=MINO(LA(I,U),LACI4N4,J)) 
LACI4+N4,J)=LACI,J) 

CONTINUE 

KK=KK+1 

IF(KK.EQ.2) GO TO 200 
NO(+4)=0 

NO(2)=0 

II=1 

DO 400 LL=1,11 

KOU(II)=1 

DO 300 I=1,4 

DO 300 J=1,4 
IF(LACI,JU).EQ LL) KOU(II)=KOUCII)+1 
CONT INUE 

IF (KOUCII).GE.7) NOCII)=LL 
IF (KOU(II).GE.7) Il=2 
CONTINUE 

RETURN 

END 


SUBROUTINE GTERP(IH,PHI,P,0,7,PG,.0G,TG,DPY ,DTY,U,UG) 
INTERPOLATES ZONAL MEAN DATA TO HEIGHT IH AND LATITUDE PHI 
DIMENSION PG(21.19),.1TG(21,19).0G(21,19) ,UG(21, 19) 


HEIGHT INDEX 
I = (IH - 15)/5 

IF (1.GT.21) I=21 

LOWER LATITUDE INDEX 

J = INT( (PHI + 100.)/10.) 

IF (J.LT.1) J = 1 

IF (J.GT.18) JU = 18 

UPPER LATITUDE INDEX 

JP = J + 1 

CHECK FOR DENSITY OR TEMPERATURE LEQ O 

CHK = DG(I,U) * TG(I.J) * DGCI,UP) * TGC(I,uUP) 
IF (CHK) 10,10,20 

P PG(1I,J) 

D DG(I.J) 

T TG(I,J) 

GO TO 30 

LATITUDE DEVIATION FROM ZONAL MFAN POSITION 

PHIF = (PHI + 100. - 10.*J)/10. 

TL= TG(I,J) + (TGCI,UP ) - TGCI,J))*PHIF 
LATITUDE INTERPOLATION 

DL= DG(I.JU} + (OGC(I,UP ) - DG(I,JU)) *PHIF 
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GTER 
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U=UG(1I,JU)+(UG(1I,JP)-UG(I,J))*PHIF GTER 
Ri = PG(1,JU)/(DG(1I,U)*TG(I.J)) GTER 
R2 = PG(I,UP)/(0G(1I,UP)*TG(I,UP)) GTER 
INTERPOLATED GAS CONSTANT GTER 
R = R1 + (R2 - R1)*PHIF GTER 
PRESSURE COMPUTED FROM INTERPOLATED GAS CONSTANT GTER 
P =DL*R*TL GTER 
D = DL GTER 
T- S211 GTER 
OP/DY FOR GEOSTOPHIC WINDS GTER 
DOPY = (PG(I,JUP) - PG(I,JU)) * 0.5 GTER 
DT/DY FOR THERMAL WINDS GTER 
OTY = (TG(I,JP) - TG(I.U)) * 0.5 GTER 
IF (ABS(PHI)-90.) 50,40,40 GTER 
DPY = O. GTER 
DTY = O. GTER 
CONT INUE GTER 
RETURN GTER 
END GTER 
SUBROUTINE INTERW(U1,V1,21,U2,V2,22,U,V,2Z) INTW 
IF ( Z1- Z2 ) 20,10,20 INTW 
u = U1 INTW 
SETS U,V = U1,V41 IF 241 = 22 INTW 
vzvi INTW 
RETURN INTW 
A = (2-21)/(Z2-21) INTW 
U = U1 + (U2-U1) * A INTW 
V = Vi + (V2-Vi) * A INTW 
LINEAR INTERPOLATION BETWEEN U1,V1 AT HEIGHT Z1 AND U2,V2 AT INTW 

HEIGHT Z2. OUTPUT IS U,V AT HEIGHT Z INTW 
RETURN INTW 
END INTW 
SUBROUTINE INTERZ(P1,01,71,21,P2,02,1T2,72,P,D,T,2Z) INTZ 
IF (Z1 - 22) 20,10,20 INTZ 
P = Pf INTZ 
D = Dt INTZ 
SETS. P..-D;.. ToS POT Fie Tk-21' = 22 INTZ 
T 274 INTZ 
RETURN INTZ 
A = (Z - Z1)/ (22 - 21) INTZ 
T = T1 + (T2 - TI)*A INTZ 
D = Di + (D2 - D1)*A * INTZ 
P = P41 + (P2 - P1) * A INTZ 
LINEAR INTERPOLATION BETWEEN P1,01,7T1 AT HEIGHT Z1 AND P2,D2,T2 INTZ 

AT HEIGHT 22 TO OUTPUT VALUES OF P,D,T AT HEIGHT 2 INTZ 
RETURN INTZ 
END INTZ 
SUBROUTINE INTER2(P1,01,71,Z1,P2,02,172,22,P,D,T,Z) INT2 
INTERPOLATES BETWEEN P1,01,T71 AT HEIGHT Z1 AND P2,02,T2 AT INT2 

HEIGHT Z2 TO OUTPUT VALUES OF P,D,T AT HFIGHT Z INT2 
CHECKS FOR T1,D1,72,D2 PRODUCT = 0, FOR GAS CONSTANT INTERPOLATIONINT2 
CHK=T1*01*T2*D2 INT2 
IF (CHK) 10,10,5 INT2 
IF (21 - 22) 20,10,20 INT2 
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.....INTERPOLATES BETWEEN 4D ARRAYS P(1I,IH),D(1I,1H),T(1,1IH) AT GRID 
LOCATIONS LATITUDE GLAT(I) LONGITUDE GLON(I). 

CLAT,.CLON = CURRENT LATITUDE ,LONGITUDE 

IZ = HEIGHT NG 
OUTPUT = P4,D4,T4, AND DERIVATIVES DPX ,DPY,DTX,DTY 


20 


30 


10 


's) 
D 


P4 

= D1 

SETS P.D,T = P1,D1,T1 IF 21=22 
T= T1 

RETURN 


IF(P1*O01*T1*P2*D2*T2.LE.0.}GO TO 30 


IF(D2*D1.LE.0.0)GO TO 30 
A=ALOG(D2/D1)/(22-21) 

LINEAR INTERPOLATION ON LOG D 
DZ= DI*EXP(A*(Z - 21)) 
A=(Z-21)/(22-2Z1) 

LINEAR INTERPOLATION ON 7 
TZ= T1 + A*(T2-T1) 
R1=P1/(D1i*T1) 

R2=P2/(02*T2) 


LINEAR INTERPOLATION ON GAS CONSTANT R 


R=(R2-R1)*A+R1 

PRESSURE FROM PERFECT GAS LAW 
= DZ * R * TZ 

DZ 

TZ 

URN 
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END 


SUBROUTINE INTER4 (€ CLAT, CLON, IZ, 


$¢ P4,D4,T4,DPX ,.OPY,DTX,DTY) 


COMMON/ IOTEMP/IOTEM1, TOTEM2,1IUS,DD,XMJD,PHI1,PHI, 


& NSAME ,DUMMY2(55) 


P, 


D, 


qT, 


NUMBER OF 4D GRID POSITIONS 


COMMON/C4/GLAT(16),GLON( 16) ,NG, DUMMY (2499) 


COMMON/CHIC/LA(4,4),NB(2), IWSYM,USH, VSH, DUSH,DVSH 
P(16,26),D(16,26),T(16,26),LAX(16) 


DIMENSION 
IWSYM=ICHAR(’ °) 

ICHK = O 

HEIGHT INDEX = HEIGHT + 1 
Id = IZ + 1 

IF (ICHK.GT.1) GO TO 220 
IF (NG.GT.9) GO TO 100 
NG = 9 MEANS POLAR GRID 


XLON=CLON 
00 10 I=10,16,1 
P(IE,IH) = P(9,IH) 
O(1,I1H) = D(9,IH) 
TCI,IH) = T(9,I1H) 
GLAT(I) = GLAT(9) 
1=10--16 ALL AT 90 DEG 
GLON(I) = GLON(I-8) 
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LOWER RIGHT INTERPOLATION INDEX INTA 
IB = INT(CLON/45) + 1 INT4 
LOWER LEFT INTERPOLATION INDEX INT4 
IA =IBt+1 INT4 
IF (IA.GT.8) IA = IA-8 INT4 
POSITION OUTSIDE POLAR GRID INT4 
IF (ABS(CLAT).LT.75.)}) GO TO 20 INT4 
UPPER LEFT INTERPOLATION INDEX INT4 
IC = IA +8 INT4 
UPPER RIGHT INTERPOLATION INDEX INT4 
ID = IB + 8 INT4 
GO TO 300 INT4 
IF (NSAME .EQ.1) NSAME=2 INT4 
CALL GEN4D INT4 
ICHK = ICHK + 1 INT4 
GO TO 5 INT4 
XLON = CLON INT4 
DO 105 I = 1,4 INT4 
DO 105 J = 1,4 INT4 
116 = 4*(I-1) + J INT4 
LAX(I16) = LA(I,uV) INT4 
CONTINUE INT4 
IF (XLON-GLON(1).GT.180.)XLON = CLON - 360, INT4 
CHECKS FOR POSITION WITHIN 16 POINT GRID 110=GOOD. 200=POSITION INT4 
OUTSIDE GRID. INT4 

IF (CLAT.GE.GLAT(1) .AND. CLAT.LT.GLAT(16) .AND. XLON.LE.GLON(1) INT4 
.AND.XLON.GT.GLON(16)) GO TO 110 INT4 
GO TO 200 INT4 
NDOL=5 INT4 
IF(ABS(CLAT).LT.18) NOL=12 INT4 
IA = 1 + INT((GLON(1) - XLON) / 5) INT4 
IA = LOWER LEFT (REFERENCE) INTERPOLATION INDEX INT4 
IA = IA + 4 * INT((CLAT ~ GLAT(1)) / NOL) INTA 
LOWER RIGHT INTERPOLATION INDEX INT4 
IB = IA + 1 INT4 
UPPER LEFT INTERPOLATION INDEX INT4 
1c = IA + 4 INI4 
UPPER RIGHT INTERPOLATION INDEX INT4 
ID = IA + 5 INT4 
GO TO 300 INT4 
IF (NSAME .EQ. 1)NSAME=2 INT4 
CALL GEN4D INT4 
ICHK = ICHK + 1 INT4 
GO TO 5 INT4 
WRITE(6, 250) INT4 
FORMAT(’ UNABLE TO GENERATE 4-D GRID. TOG MANY RETRIES IN INTER4’ )INT4 
P4=0. INT4 
D4=0. INT4 
T4=0. INT4 
RETURN INT4 
INTERPOLATION FOR POSITION INSIDE 1G POINT GRID OR POLAR GRID INT4 
300 CALL INTLL(P,IA,IB,IC,10,P4,GLAT,GLON,CLAT,XLON, IH) INT4 
CALL INTLL(D,IA,IB,IC,10,D4,GLAT,GLON,CLAT, XLON, IH) INT4 
CALL INTLL(T,IA,IB,I1C,10,T4,GLAT,GLON,CLAT,XLON, IH) INT4 
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RELATIVE LONGITUDE DISPLACEMENT FROM REFERENCE POSITION (IA) 
= (XLON - GLON(IA))/(GLON(IB) - GLON(IA)) 

RELATIVE LATITUDE DISPLACEMENT FROM REFERENCE POSITION(IA) 
(CLAT - GLAT(IA))/(GLATCIC) - GLAT(IA)) 


DLON 


DLAT 


DPX=P(IB,IH)-P(IA,IH) 


DP/DX FOR GEOSTROPHIC WIND EQUATIONS 

DPX = DPX + (P(ID,IH) - PCIC,IH) - DPX)*DLAT 
DTX = TCIB,IH) - TCIA,IH) 

DT/DX FOR THERMAL WIND EQUATIONS 

OTX = DIX + (TCID.IH) - TEIC,IH) - OTX)*DLAT 
DPY = P(IC,IH) - PCIA,IH) 

DP/DY FOR GEOSTROPHIC WIND EQUATIONS 

OPY = DPY + (P(ID,IH) - P(IB,IH) - OPY)*DLON 
DTY = T(IC,IH) - TCIA,IH) 

DT/DY FOR THERMAL WIND EQUATIONS 

DTY = DOTY + (TCID,IH) - TCIB,IH) - DTY)*DLON 
IF(NG.GT.9) GO TO 315 

DPX=DPX/9. 

DTX=DTX/9. 


DPY=SIGN(OPY/3. ,CLAT) 

OTY=SIGN(OTY/3.,CLAT) 

IF(ABS(CLAT).GT.18) GO TO 312 

DPY=DPY*5./12 

DTY=DTY*5./12 

RETURN 

END 

SUBROUTINE INTLL(F.IA,IB,IC,ID,FLL,GLAT,GLON,CLAT,CLON, IH) 

INTERPOLATES FUNCTION (ARRAY) F FROM VALUES OF GLAT AND GLON AT 
INDEX VALUES IA, IB, IC, ID TO OUTPUT VALUE FLL AT HEIGHT IH 
AND POSITION CLAT, CLON 

DIMENSION F(16,26),GLAT(16),GLON(16) 


-.NORMALIZES LONGITUDE DISPLACEMENT 


IF(FC(IA,IH)*F(IB,IH)*F (IC, IH)*F(ID,IH)) 20, 10,20 

FLL=O. 

RETURN 

X=(CLON-GLON(IB))/(GLON(IA)-GLON( IB) ) 

NORMALIZES LATITUDE DISPLACEMENT 

Y=(CLAT-GLAT(IA))/(GLAT(CIC)-GLAT(IA)) 

TWO DIMENSIONAL INTERPOLATION 

FLL=FCIB.IH)+(F CID, IH)-FCIB,IH))*Y+C(F CIA, IH)-FCIB,IH))*xX 
+(FCIC.IH)-FCIA,IH})-FCID, IH)+F (1B, IH) )*X*Y 

RETURN 

END 

SUBROUTINE INTRP4 (LALON) 


SUBROUTINE TO INTERPOLATE VALUES 
DIMENSION XLL(4).YLL(4).XC(4),YC(4) 
COMMON/ INT/D(208,5),1G(5),0xX¥(2),DLA(4),DLO(4) 
DEGRAD=3.14159/180. 


LALO=IABS(LALON) 
L1=LALO/ 10000 
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4 (D(I,4)-D(1,3)-D(1,2)+D(1,1) 


L2=LALO-L1* 10000 
XL=L1/10. 

YL=L2/10. 

IF (I1G(5)-2) 30,20,10 
IF (1G6(5)-3) 30,30,50 


INTERPOLATE FROM NMC GRID 


CONT INUE 

DO 25 L=1,26 

DO 22 J=1,4 

IF (O(tL,J).LT.0.01) GO TO 25 

DO 24 K=1,8 

I=(K-1)*264L 
D(I,5)=(1.-Dx¥(2))*(€1.-OxXY¥(1))*DC1L,1)4+DxXY¥(1)*O(1,2)) 


1 +Dxvy(2)*(C€C4.-DXY¥(4))*001,3))4+DxXY(1)*0(1,4)) 


CONTINUE 
CONT INUE 
RETURN 


INTERPOLATE FROM EQUATION FOR SOUTHERN HEMISPHERE GRID 


CONT INUE 

DO 32 J=1,2 

XLL(J)=DLA(U) 

YLL(J)=DLO(U) 

IF ((YL.GE.355.).AND.(YLL(JU).LT.0.01)) YLL(U)=360. 

CONTINUE 

X=(YLL(1)-YL)/5. 

Y=(XL-XLL(1))/5. 

IF (1G(5).£Q.3) Y=-Y 

DO 38 L=1,26 

DO 36 J=1,4 

IF (D(L,J).LT.O.O01t) GO TO 38 

DO 37 K=1,8 

I=(K-14)*264+L 

D(I,5)=D(I,1)4+X*(D(1,2)-D(1,1))+¥*(D01,3)-D(1,1))+x*Y¥* 
I ) 

CONTINUE 

CONT INUE 

RETURN 


INTERPOLATE FROM ACROSS GRIDS 


CONTINUE 

IF (1G(5).NE.1133) GO TO 55 
1G(5)=3 

GO TO 30 

CONT INUE 

IF (IG(5).NE.333) GO TO 60 
DLO(1)=(DLO(2)+DIL0(3))/2. 
DO 52 1=1,208 

D(I,4)=D(1,3) 

DLA(4)=DLA(3) 
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DLO(4)=DLO(3) 
CONTINUE 

DO 62 I=1,4 
XLL(I)=DLA(T) 
YLL(1)=DLOCI) 


IF ((YL.GT.350.).AND.(YLL(I). LT.0.01)) YLL(I)=360. 


CONTINUE 

ITH=O0 

X=YLLOC1)-Y¥L 

Y=XL-XLL(1) 

CONTINUE 

DO 65 I=2,4 

XC(I)=YLLC1)-YLLCT) 
YC(1)=XLL(I)-XLL(1) 

TH2=3.14159/4 

TH3=3.14159/4 

IF (ABS(XC(2)).GT.0.01) TH2=ATAN(YC(2)/XC(2)) 
IF (ABSCYC(3)).GT.0.01) TH3=ATAN(XC(3)/YC(3)) 
IF (XC(2).LT7.0.) TH2=3.14159+TH2 

IF (xcC(3).LT.0.) TH3=3.141594+TH3 
DNN=COS(TH2+TH3) 

IF (ABS(DNN).GT.0.001) GO TO 66 
ITH=ITH+1 

IF (ITH.EQ.2) GO TO 66 

XLLOE3)=XLL(4) 

YLL(3)=Y¥LL(4) 

O00 61 I=1,208 

O(1,3)=D(1I,4) 

GO TO 63 

CONTINUE 

ZA=SORT(XC(2)**2+YC(2)**2) 

IF (ITH.LT.2) GO TO 69 
Z=SQRT(X**2+4+Y**2) 

E=0. 

Z24=0. 

GO TO 71 

CONT INUE 

EB=SQRT(XC(3)**2+YC(3)**2) 
24=(XC(4)*COS(TH3)-YC(4)*SIN(TH3))/DNN 
E4=(¥C(4)*COS(TH2)-xC(4)*SIN(TH2))/DONN 
Z2=(X*COS(TH3)-Y*SIN(TH3) )/ONN 
E=(Y*COS(TH2)-X*SIN(TH2))/DONN 

B=0. 

C=O. 

DD=-oO. 


CONTINUE 
DO 70 L=1,26 

DO 68 J=1,4 

IF (D(L,J).tLT.0.01) GO TO 70 

DO 67 K=1,8 

T=(K-1)*264L 

A=D(I, 1) 

IF (2A.GT.0.01) B=(O0(1,2)-D(1,1))/ZA 
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IF (€B.Gr.0.01) Cz=(N€1,9)-DC(1I,1))/EB 


IF (CABS(71).6T.0.01). AND. (ABS(E4).GT.0.901)) 


t DD=(NC1,4)-A-BtZ4-Ct+E4)/(Z47E4) 

Ol AARC EDO ee 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE TNTRUV(UR,VR,H, PHI ,SUH,SVH) 


PHI (DEGREES). FROM UR AND VR ARRAYS 
DIMENSION UR(29,19),VR(29, 19) 


be pcces I - LOWER tHETGHT INDEX 


hagas PHI1 - LOWER LATITUDE FOR UR AND VR ARRAY VALUES 


T= 1+ TNT (t1)/5 

IF (H.GE 125) 1=25'CINTOH)-120)/20 
If (1.Q1T. 29) 129 
UPPER HEIGHT INDEX 
IP=I+1 

IF (IP.GT.29) IP=29 
LOWER LATITUDE INDEX 
J=INT(PHI+100. )/10 
UPPER LATITUDE INDEX 
JP=J+1 

IF (uP.GT.19) JP=19 


PHI 1=-100.4+10.*d 


PHI2=-100.+10. *uUP 
IF (1,GT.25) GO TO 10 


LOWER HEIGHT FOR UR AND VR ARRAY VALUES 


Z1=5.*(1I-1) 

GO TO 20 
Z1=20.*(1-19) 
IF(IP.GT.25) GO TO 30 


UPPER HEIGHT FOR UR AND VR ARRAY VALUES 


72=5.*(IP-1) 

GO TO 40 

72=20.*(IP-19) 

INTERPOLATE ON LATITUDE AT LOWER HEIGHT 


40 CALL INTERW(UR(I,JU),VR(1,J0),PHI1,UR(I,JuP),VR(I,JP),PHI2,U1.V1, 


FINDS RANDOM WIND STANDARD DEVIATION AT HEIGHT t! (KM), LATITUDE 


PHI2 - UPPER LATITUDE FOR UR AND VR ARRAY VALUES 


$PHI ) 


CALL INTERW(URCIP,. J), VRCIP, J), PHI. URCIP,JP),VRCIP,UP),PHI2,U2,V2. 


INTERPOLATE ON LATITUDE AT UPPER HEIGHT 


$PHI ) 
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= 2 © # © 


INTERPOLATE ON HEIGHT 


CALL INTERW(U1,V1,21,U2,V2,Z2, SUH, SVH,H) 
RETURN 

END 

SUBROUTINE INTRW(WR,H, SWH) 


FINDS RANDOM WIND STANDARD DEVIATION AT HEIGHT H (KM), 


FROM WR ARRAYS 


DIMENSION WR(29) 


I 
I 
I 


- LOWER HEIGHT INDEX 
=4+INT(H)/5 
F (H.GE.125) 1225+(INT(H)-120)/20 


INTFR 
INTR 
INTR 
INTR 
INTR 
INTR 
INTR 
INTR 
INUV 
INUV 
INUV 
INUYV 
INUV 
TNUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUYV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
INUV 
TNUV 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 


120 
121 
122 
123 
124 
125 


—- ane 
ON 
“oO 


OmrN OA TS&wWQDh — 


29 


30 


10 
20 


le) 


40 


30 


40 CALL INTERW(UR(I,J).VRCT,J),PRHI1,UR(T,JP),VR(I,JP),PHI2,U1,V1, 


IF (1.GT.29) 1=29 
UPPER HEIGHT INDEX 
IP=I+1 
IF (IP.GT.29)} IP=29 
IF (I1.EQ.IP) GO TO 50 
IF (I.GT.25) GO TO 10 
LOWER HEIGHT FOR WR ARRAY VALUES 
Z71=5.*(1-1) 
GO TO 20 
Z1=20.*(1-19) 
IF(IP.GT.25) GO TO 30 
UPPER HEIGHT FOR WR ARRAY VALUES 
Z2=5.*(IP-1) 
GO TO 40 
72=20.*(1IP-19) 
INTERPOLATE ON HEIGHT 
A=(H-21)/(22-2Z1) 
SWH=WR(I) + (WRCIP)-WR(I))*A 
RETURN 
SWH=WR(I) 
RETURN 
END 
SUBROUTINE INTR25(UR,VR,H,PHI , SUH, SVH) 


FINDS LARGE SCALE FRACTIONAL VARIANCES AT HEIGHT H (KM),LATITUDE 


PHI (DEGREES), FROM UR AND VR ARRAYS 

DIMENSION UR(25, 10),VR(25, 10) 

I - LOWER HEIGHT INDEX 

IF (H.LT.95.) I = t + INT(H) / 5 

IF (H.GE.95.) 1=19+(CINT(H)-80)/20 

IF (1.GT.25) I = 25 
UPPER HEIGHT INDEX 

IP=I+1 

IF (IP .GT.25) IP=25 

LOWER LATITUDE INDEX 

J=INT(PHI+110.)/20 

UPPER LATITUDE INDEX 

JP=J+ 1 

IF (JP.GT.10) JP=10 

PHI4 - LOWER LATITUDE FOR UR AND VR ARRAY VALUES 
PHI 1=-110.4+20.*J 

PHI2 - UPPER LATITUDE FOR UR AND VR ARRAY VALUES 
PHI2=-110.4+20.*uUP 

IF (L.GT.198) GO TO 10 

LOWER HEIGHT FOR UR AND VR ARRAY VALUES 
21=5.*(1-1) 

GO TO 20 

Z1=20.*(1-15) 

IF (IP.GT.18) GO TO 30 

UPPER HEIGHT FOR UR AND VR ARRAY VALUES 
Z2=5.*(IP-1) 

GO TO 40 

Z2 = 20. * (IP - 15) 

INTERPOLATE ON LATITUDE AT LOWER HEIGHT 


INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
INRW 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25 
IN25S 
IN25 
IN25 
IN25 
IN25 


$PHI ) 


INTERPOLATE ON LATITUDE AT UPPER HEIGHT 
CALL INTERW(UR(IP.J), VR(IP.J).PHI1,URCIP, JP). VRCIP,JP),PHI2,U2,V2,1N25 


$PHI ) 


INTERPOLATE ON HEJGHT 
CALL INTERW(U1,V1,21,.U2,V2,22,SUH, SVH,H) 


RETURN 
END 


SUBROUTINE JAC(Z.T, TZ, AN, A02,A0,AA,AHE ,AH,EM,OENS,DL) 


IN25 
IN25 


IN25 
IN25 
IN25 
IN25 
IN25 
JAC 

JAC 


Ce SS EI HIE TTT EH TT HEEL REEL ELSE SSAC 


Ct 
Cc ** 
C** 
C** 
C** 
C* * 
Cr* 
C** 


‘JAC’ CALCULATES THE TEMPERATURE TZ , THE TOTAL DENSITY DENS AND 
ITS LOGARITHM DOL, THE MEAN MOLECULAR WEIGHT EM, THE INDIVIDUAL 
SPECIE NUMBER DENSITIES FOR N, O02, O, A, HE AND H (EACH PRECEDED 
WITH AN ’A’) AT ALTITUDE Z GIVEN THE EXOSPHERIC TEMPERATURE ¥: 
THIS SUBROUTINE USES THE SUBROUTINE ‘GAUSS’ AND THE FUNCTION 


SUBPROGRAMS ‘TEMP’ AND ‘MOLWT’. 


** JAC 
** JAC 
** JAC 
** JAC 
YF JAC 
¥t JAC 
t*¥ JAC 
** JAC 


CORR III IOI IOI RIOT TERETE TT EEE TER EE LT ERET ELEN TS *+** JAC 


C** 
C** 


C 


PARAMETER (AV = 6.02257E23) 
PARAMETER (QN = .78110) 
PARAMETER (QO2 = .20955) 
PARAMETER (QA = ,009343) 
PARAMETER (QHE = 1.289E-05) 
PARAMETER (RGAS = 8.31432) 
PARAMETER (PI = 3.14159265) 
PARAMETER (TO = 183.) 


DIMENSION ALPHA(G) , EI(6) , DI(6) , DIT(6) 


REAL 


DATA 
DATA 


TX = 444.3807 + .02385 * T - 392.8292 * ExP ( -.0021357 * T ) 
A2 = 2. * (T-TX) / PI 
TXTO = TX - TO 
Ti = 4.9 * TXTO /* 35. 
T3 = -4.7 * TXTO / ( 35.**3 ) 
T4 = -0.8 * TXTO / ( 35.**4 ) 
Tz = TEMP ( 2. TX , T1,73 =, 74, A2 ) 
SECTION 1 
D = AMINi ( Z , 105. ) 
INTEGRATE GM/T FROM 90 TO MINIMUM OF Z OR 105 KM :- 


MOL WT 
ALPHA / 0.0 , 0.0, 0.90 , 0.0 , ~.380 , 0.0 / 
£1/28.0134,31.9988, 15.9994,39.948,4.9026,1.00797/ 


JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 


31 


32 


CALL GA 


C THE NUMBER 

C AT 90 KM. 
EM 
TO 


DENS 


FACTO 
PAR = 
FACTO 


USS(D.41,R,TX,.T1.73,14,A2) 
2.1926E-8 = DENSITY X TEMPERATURE/MEAN MOLECULAR WEIGHT 


MOLWT ( D ) 
TEMP ( O , TX , T1 ,739 , T4, AQ ) 


= 2.1926E-8 * EM * EXP( -R / RGAS ) / TD 
R = AV * DENS 


FACTOR / EM 
R = FACTOR / 28.96 


C FOR ALTITUDES BELOW AND AT 105 KM CALCULATE THE INDIVIDUAL SPECIE 


C NUMBER DEN 
IF { 


DL = 
AN 
AA 
AHE 
AO 
AO2 
AH 


| | 


C 
C** RETURN TO 
C 


SITIES FROM THE MEAN MOLECULAR WEIGHT AND TOTAL DENSITY. 
Z. LE. 105 ) THEN 


ALOG1O ( DENS ) 

ALOG1O ( QN * FACTOR ) 

ALOG1O ( QA * FACTOR ) 

ALOG1O ( QHE * FACTOR ) 

ALOG1O ( 2. * PAR * ( 1.-EM / 28.96 ) ) 

ALOG10 ( PAR * ( EM * ( 4.+#Q02 ) / 28.96-1. ) ) 
O. 


CALLING PROGRAM 


RETURN 


END I 


F 


C** SECTION 2 : THIS SECTION IS ONLY PERFORMED FOR ALTITUDES 


aaa ABOVE 105 KM. 


C NOTE THAT HAVING REACHED THIS SECTION MEANS THAT D IN SECTION 1 IS 


C 105 KM. 


C CALCULATE INDIVIDUAL SPECIE NUMBER DENSITIES FROM THE TOTAL OENSITY 


C AND MEAN MOLECULAR WEIGHT AT 105 KM ALTITUDE. 


C INTEGRATE 


DI( +) QN * FACTOR 


DI(2) = PAR * (EM * (1.4902) / 28.96-1.) 
DI(3) = 2. * PAR * (1.- EM / 28.96) 
DI(4) = QA * FACTOR 

DI(5) = QHE * FACTOR 


G/T FROM 105 KM TO 2 KM :- 


JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 


JAC 


JAC 
JAC 
JAC 
JAC 


CALL GAUSS(Z,2,R,1X,71.T3,74,A2) 
DO. 41° 2 sf 3-5 
DIT(I) = OI1(1)*(TD/TZ)**(14.4+ALPHACT))*EXP(-EI(1)*R/RGAS) 


IF ( DIT(I). LE. O. ) DIT(T) = 1.E-6 
41 CONT INUE 


C** THIS SECTION CALCULATES ATOMIC HYDROGEN DENSITIES ABOVE SOO KM 
C** ALTITUDE. BELOW THIS ALTITUDE . H DENSITIES ARE SET TO 10**-6., 


crs SECTION 3 
Cee =e -- 


IF ( Z .GT. 500. ) THEN 


Ai = 500. 
S:.=: FEMP “CAD 6 TK OT Gs 14 4 TS a AZ DD 


DI(6)=10.**(73.13-39.4*AL0G10(S)45.5*ALOG10(S)*ALOG10(S)) 
CALL GAUSS(Z,7,R,1TX,T1,73,7T4,A2) 
DIT(6) = DI(6) * (S/TZ) * EXP ( -EI(6) * R / RGAS ) 
ELSE 
DIT (6) = 1.0 


END IF 


C FOR ALTITUDES GREATER THAN 105 KM , CALCULATE TOTAL DENSITY AND 
C MEAN MOLECULAR WEIGHT FROM INDIVIDUAL SPECIE NUMBER DENSITIES. 


DENS=O 
DO: 42 “T= 1 Ss 
DENS = DENS + EI(I) * DIT(1I) / AV 
42 CONTINUE 


EM=DENS*AV/(DIT(1)4DIT(2)4+DIT(3)4D0IT(4)+DIT(5)+DIT(6)) 


DL = ALOGiO (DENS) 


AN = ALOGIO(DIT(1)) 
AO2 = ALOGIO(DIT(2)) 
AO = ALOGIO(DIT(3)) 
AA = ALOGIO(DIT(4)) 
AHE = ALOGIO(DIT(5)) 
AH = ALOGIO(DIT(6)) 


JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 
JAC 


101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
1114 
172 
113 
114 
115 
116 
117 
118 
119 
120 
124 
122 
123 
124 
125 
126 
127 
128 
A239 
130 
131 
132 
133 
134 
135 
136 
137 
1398 
{99 
140 
141 
142 
143 
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 


33 


34 


RETURN 


END 


SUBROUTINE JACCH(Z,PHIR, THET,PH,DH, TH) 
COMMON/TOTEMP/101,102,103,05,06,07,D8, 
& 109,010,011,012,D13,D14,015,016,017,D18,D19, 


& M, 


& THR,MIN,028 ,0D29,030,D31,D32,DUMMY(24) 
REAL INDATA( 12) ,OUTDATA( 12) 
PI=3.14159265 


INDATA(1)=2Z 


INDATA(2)=PHIR* 180. /PI 
INDATA(3)=-THET 
INDATA(4)=FLOAT(IYR) 
INDATA(5)=FLOAT(M) 
INDATA(6)=FLOAT( IDA) 
INDATA(7)=FLOAT(C IHR) 
INDATA(8)=FLOAT (MIN) 
INDATA(9)=2. 
INDATA(10)=F tO 

INDATA( 11)=F 108 
INDATA( t2)=AP 

CALL J7O(CINDATA,GUTDATA) 
PH=OUTDATA( 12) 
DH=OQUTDATA( 10) 
TH=OUTDATA(2) 


RETURN 
END 


SUBROUTINE J70 ( INDATA, 


OUT 


DATA ) 


IDA, fTYR,020,021,022,D23,D24,0D25,D26 ,027,F10,F10B, AP, 


JAC 

JAC 

JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
JACC 
J70 

J7O 


CERT E RRR ERE EEE EE ETH EE EEE EEA E REET EEE REET RE REECE EEEH EERE RENE HERE KE HE OITO 


Ct 
C*+ 
Cer 
Ce 
Cx 
O55 
C** 
C+ 
Cr* 
Ce 
C** 
C** 
c+ 
C** 
C++ 
C** 
Cr+ 
C++ 
C*r* 
Ce* 
Ctr 
C* 
Cc * * 
C+* 


11 
F100 


F108 


GI 


J7O DEVELOPED FROM J7OMM BY 


MIKE P. HICKEY 
UNIVERSITIES SPACE RESEARCH ASSOCIATION 
AT 
NASA / MARSHALL SPACE FLIGHT CENTER, £044, 
HUNTSVILLE, ALABAMA, 35812, USA. 
TEL. (205) 544-5692 
INPUT DATA: 
~ ALTITUDE = INDATA (1 
- LATITUDE = INDATA (2 
- LONGITUDE = INDATA (3 
YEAR (YY) = INDATA (4 
MONTH (MM) = INDATA (5 
DAY (DD) = INDATA (6 
-~- HOUR (HH) = INDATA (7 
MINS (MM) = INDATA (8 
GEOMAGNETIC INDEX = INDATA (9 
SOLAR RADIO NOISE FLUX = INDATA (1 
- 162-DAY AVERAGE FIO = INDATA (1 
GEOMAGNETIC ACTIVITY INDEX = INDATA (1 


RD me ee ee ee ee ee ee ee 


a a 


**J70 
*+*J70 
** J70O 
**JI70 
**J70 
*#*J70 
**JI70 
**JU70 
**J'7O 
**J70 
**J70 
** 570 
**J70 
** 70 
** N70 
** 570 
** ITO 
**J70 
**I7O 
*I7O 
**J70 
** ITO 
= eTO 
**UTO 


155 
156 


C** 
C** 
C=* 
C*r 
C** 
Cs 
C** 
C** 
ct. 
C** 
Cex 
C** 
C** 
C** 
C** 
Cee 
C* * 
C** 
Cc * * 


NB. 


DENS-- 


FM-- AVERAGE MOLECULAR WEIGHT 
TOTAL DENSITY 
DL-- LOG1O ( TOTAL DENSITY ) 


A(1)-- N2 NUMBER 
A(2)-- O02 NUMBER 
A(3)-- QO NUMBER 
A(4)-- A NUMBER 
A(5)-- HE NUMBER 
A(6)-- H NUMBER 


OUTPUT 


—— ee eee 


T-- EXOSPHERIC TEMPERATURE 
TZZ-~ TEMPERATURE AT ALTITUDE 2 


DENSITY 
DENSITY 
DENSITY 
DENSITY 
DENSITY 
DENSITY 


P-- TOTAL PRESSURE 


INPUT THROUGH! ARRAY 
OUTPUT THROUGH ARRAY 


DATA: 


Yo oto it ou& 


1 ee | ee) 


lt oH 


‘INDATA’ 
‘OUTDATA’ 


OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 
OUTDATA 


OUTDATA { 
QUTDATA ( 
QUTDATA ( 


fa fin pen pean el pi A em, 


**U70 
*+*J70 
** STO 
¥ *JI70 
+*#J70 
**J70 
**J70 
**J70 
Tea 7TO 
*¥*JI70 
** ITO 
** 570 
**J70 
**I7O 
at 36, 
£10 
** ITO 
eras TO 
*+J70 


CR OIF IO OI III OIG TOO IO IT IIH EH HATER EE YY ITO 


C CALCULATIONS PERFORMED FOR ONLY ONE LATITUDE 


PARAMETER (RGAS = 8.31432E3) 


PARAMETER (BFH 


DIMENSION A 


REAL 


INDATA 


C AND ONE ALTITUDE 


C 


C** SET PARAMETERS TO 


C 


Z 
XLAT 
XLNG 
IYR 
MN 
IDA 
IHR 
MIN 
I 1 
F10 
F 10B 
GI 


CALL 


CALL 


CALL 


INDATA 
INDATA 
INDATA 
INT ( 

INT { 
INT ( 
INT { 
INT ( 
INT { 

INDATA 
INDATA 
INDATA 


nou uw HW UW & ht nH ot ot on a 


JAC(Z,TE.TZ.A( 1), A(2),A03),A04),A05),A06),EM,DENS,DL) 


DENLG 


( 6 


C2?) 


440.0) 
) 


QOUTDATA ( 


INDATA VALUES 


(1) 
(2) 
(3) 
INDATA 
INDATA 
INDATA 
INDATA 
INDATA 
INDATA 
(10) 
(11) 
(12) 


TINE ( F100 , F108 


= OQ. 


+ 


. GI, 


, ONE LONGITUDE 


TME(MN.IDA,IYR, IHR,MIN, XLAT,XLNG,SDA,SHA,DD,DY ) 


SDA , SHA 


J7O 
J70 
J7O 
J70 
J7O 
J70 
J7O0 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J7O 
J70 
J70 
J7O 
J70 
J7O0 
J7O 
J7O 
) J7O 
J7O 
J70 
J7O 
J7O 


36 


C 
cry 
C 


C 


80 


DUMMY 
DEN 


DL. 
DL 


IF ( Z .LE. 170. ) THEN 
CALL SLV ( DUMMY . Z , XLAT , DD ) 
DENLG = DUMMY 

END IF 


‘FAIR’ HELIUM NUMBER DENSITY BETWEEN BASE FAIRING HEIGHT ( BFH ) 
AND 500 KM. 


IF ( 2. GE. 500. ) THEN 
CALL SLVH ( DEN , ACS) , XLAT , SDA ) 
DL = DEN 
ELSE IF ( 2 .GT. BFH ) TH 
DHEL{1 = A 
DHEL2 = A 
DLG1 = DL 
DLG2 = DL 
CALL SLVH ( DLG2 
IH = 2 
CALL FAIRS(DHEL1,OHEL2,DLG1,OLG2,IH,FDHEL.FDLG) 
DL = FDOLG 
A (5 ) = FOHEL 


, DHEL2 , XLAT , SDA ) 


END IF 


BDL = DL + DENLG 
DENS = #t0.7*DL 
XLAT = XLAT * 57.29577951 


FILL OUTDATA ARRAY 
OUTDATA (1) = TE 
OUTDATA (2) = TZ 


DO 8O I = it, & 
OUTDATA (1+2) = 41.E6 * ( 10. ** A(T) ) 
CONT INUE 


OUTDATA (9) = EM 

OUTDATA (10) = DENS * 1000. 
OUTDATA (11) = DL 

P = OUTDATA (10) * RGAS * T2 / EM 
OUTDATA (12) = P 


RETURN 
END 
REAL FUNCTION MOLWT( A ) 


CREEK ERR KR ORK EO RE Re RO KOK KR RR Ee ERO ROKR 


C* * 
Cer 
C** 


“MOLWT’ CALCULATES THE MOLECULAR WEIGHT FOR ALTITUDES BETWEEN 
30 AND 105 KM ACCORDING TO EQUATION (1) OF SAO REPORT 313. 


> 
“+ 
* * 


J70 
J70O 
J70 
J70 
J7O 
J70 
J7O0 
J70 
J7O 
JTO 
J70 
J70O 
J7O 
J70 
J7O 
J7O 
J70 
J7O 
J7O 
J7O 
J70 
J7O0 
J70 
J7O 
J7O 
J70 
J70O 
J70 
J7O 
J7O 
J70 
J70 
J7O 
J7O 
J70O 
J7O 
J7O 
J70O 
J7O 
J7O 
J7O 
J70O 
J7O 
J70 
J7O0 
J70 
J70 
J7O 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 


112 
113 
114 
115 
116 
VAT 
118 
119 
120 
12 1 
122 
123 
124 
125 
126 
127 
128 


On hr WA -- 


C** OTHERWISE, MOLWT IS SET TO UNITY. a 


cC*« 


** 


Ce KE TNH I ER BHT TEETER ERR EAH ETT ETT ES 


& 


& 


DIMENSION 8 (7) 


DATA B / 28.15204 , -0.085586, 1.284E-4, -1.0056E-5, -1.021E-5, 
1.5044E-6, 9.9826E-8 / 
IF ( A. GT. 105. ) THEN 
MOLWT = 1. 
ELSE 
Uo = A - 100. 
MOLWT = B (1) 
HO: 4. be 
MOLWT = MOLWT + 8 (I) * U ** ( I-t ) 
CONTINUE 
END IF 
ENO 


SUBROUTINE PDTUV (PSP,DSP,TSP,USP,VSP,CLAT,CLON,IH,PS,DS,TS, 

DPX ,.DPY ,DTX ,OTY,US,VS) 

INTERPOLATES STATIONARY PERTURBATIONS ON LATITUDE AND LONGITUDE 
AT HEIGHT IH 

DIMENSION PSP( 15.19, 18).DSP( 15,19. 18). TSP( 15,19, 18),USP(15,19, 18), 
VSP( 15,19. 18) 

HEIGHT INDEX K 

K=(IH-15)/5 

XLON = CLON 

LOWER LONGITUDE INDEX J 


MOWT 
MOWT 
MOW! 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
MOWT 
PDTU 
PDOTU 
PDTU 
POTU 
POTU 
PDTU 
PDTU 
PDTU 
POTU 
PDTU 


DLON - RELATIVE LONGITUDE DEVIATION FROM CORNER REFERENCE LOCATIONPDTU 


IF(XLON .LE. 180. )THEN 

JZ INT( (200. -XLON)/20. ) 
DLON=(200. -XLON-20.%*J)/20. 
ELSE 

J=INT((560.-XLON)/20. ) 
DLON=(560.-XLON-20. *J)/20. 

END IF 
UPPER LONGITUDE INDEX uP 

JP = J+ 

IF(UP .GT. 18) JP=1 

LOWER LATITUDE INDEX I 
I=INT((CLAT+100. )/10.) 

UPPER LATITUDE INDEX IP 

IP = I+1 

IF(IP .GT. 19)1P=19 

DLAT - RELATIVE LATITUDE DEVIATION FROM CORNER REFERENCE LOCATION 
DLAT=(CLAT -10.*1+100. )/10. 
PRESSURE LAT-LON INTERPOLATION 


PDTU 
POTU 
POTUL 
POTU 
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POTU 
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PDTU 


37 


38 


oe 4 + 


. 8 6 


~“ 2 © & 


PS=PSP(K,I,JU)+(PSP(K,IP,J)-PSP(K,I,U))*DLAT+(PSP(K,I,JUP)-PSP(K,1I.JUPOTU 


1))*DLON+(PSP(K,IP.UP)-PSP(K,I.JP)-PSP(K.IP,J)+PSP(K,1,J))*DLAT* 


2DLON 
DENS 


ITY LAT-LON INTERPOLATION 


POTU 
POTU 
PDTU 


DS=DSP(K,1,JU)+(OSP(K,IP,J)-DSP(K,1I,JU))*DLAT+(DSP(K,I,JUP)-OSP(K,I,UPDTU 


1))*DLON+(DSP(K,IP,JP)-DSP(K,1I,.JUP)-OSP(K,IP,J)+O0SP(K,I,JU))*DLAT* 


2DLON 
TEMP 


ERATURE LAT-LON INTERPOLATION 


PDTU 
PDTU 
PDTU 


TS=TSP(K,1,J)+(TSP(K,IP,JU)-TSP(K,I,0))*DOLAT+(TSP(K,I,UP)-TSP(K.1I.JPOTU 


1))*DLON+(TSP(K,IP,JP)-TSP(K,I,UP)-TSP(K.IP,J)+TSP(K,I,J))*DLAT* 


2DLON 
ZONA 


L WIND LAT-LON INTERPOLATION 


POTU 
POTU 
POTU 


US=USP(K,I,JU}+(USP(K,IP,J)-USP(K,1.U))*DLAT+(USP(K,1I,JUP)-USP(K,1I,JUPOTU 


1))*DLON+(USP(K,IP,JP)-USP(K,1L,UP)-USP(K, IP, J)+USP(K,I,U))*DLAT* 


20t0N 


MERIDIONAL WIND LAT-LON INTERPOLATION 
VS=VSP(K,1I,JU)+(VSP(K,IP,J)-VSP(K,1,JU))*DLAT+(VSP(K,I,UP)-VSP(K,I,JUPDTU 


1))*DLON+ (VSP(K,IP,JP)-VSP(K.I,.JUP)-VSP(K,IP,JU)+VSP(K.I,JU))*OLAT* 


2DLON 
..OPX - DP/DX FOR GEOSTROPHIC WINDS 

DPX = (PSP(K,1I,JU)-PSP(K,1,JUP))/4. 

DPX = DPX + ((PSP(K,IP,J)-PSP(K,IP,JP))/4. - DPX)*DLAT 
..DPY - DP/DY FOR GEOSTROPHIC WINDS 

DPY=(PSP(K,IP,J)-PSP(K,1I,JU))/2. 

DPY= DPY + ((PSP(K,IP,JP)-PSP(K,I1,JUP))/2. - DPY)*DLON 
..DTX - DT/DX FOR THERMAL WINDS 

DTX = (TSP(K,I,J)-TSP(K,1,UP))/4. 

OTX = DTX + ((CTSP(K,IP,JU)-TSP(K,IP,JUP))}/4. - DIX)*OLAT 
..DTY - OT/DY FOR THERMAL WINDS 

DTY=(TSP(K,IP,J)-TSP(K,I,J))/2. 

DTY= DTY + ((TSP(K,IP,JP)-TSP(K,I,UP))/2. - DTY)*DLON 

RETURN 

END 


SUBROUTINE PERTRB 
COMMON /IOTEMP/IOTEM1, IOTEM2,IUS,.DD,XMJD,PHI1,PHI,NSAME, 
$PL1,D0L1,TL1,SPL1,SD0L1,STL1,UL1,VL1,SUL1,SVL1,MN,IDA,IYR, 


1PH,P 


* PLON,G,R,CH,CLAT,CLON,F10,F10B,AP,IHR,MIN,NMORE ,OX,HL,VL,OZ, 


LAT, 


2B,EPS,1IOPP,LOOK,IET,FLAT,PS1.0S1,7TS1,US1,VS1,5PS1,SDS1, 


3STS1 


.5US1,SVS$1,U0S1,V0S1,UDL1,VDOL1,UDS2,VDS2,UDL2,VOL2 


COMMON /COMPER/SP2,S02,ST2,P2,02,T2,U2,V2,SU2,SV2,CP, 
1PS2,082,TS2,US2,.VS2, 
2PL2,0L2,TL2,UL2,VL2, 


3SPS$2 
4SPL2 


»9DS2,STS2,SUS2,SVS2, 
,90L2,STL2,SUL2,SVL2 


COMMON /VERT/W1,SW1,W2,SW2,WR(29) 
COMMON /WINCOM/DUM(8),T,DUMMY2(9) 
COMMON /RANOD/IX,1Y,12 


DLON 
PI = 


= ABS(CLON-PLON) 
3.1415927 


IF(DLON.GT.PI) DLON = 2.*PI - DLON 


DX 
IF ( 


= R*SORT((CLAT-PLAT)**2 + (COS(CLAT)*(DLON 
ABS(DX).LT.O.01.AND.ABS(D02Z).LT.0.01)DX = 0.01 


))**2) 
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DX IS HORIZONTAL DISTANCE BETWEEN POSITIONS PLAT.PLON AND CLAT,CLOPERT 


10 
20 


30 
40 


50 
60 


790 
80 


90 
100 


110 
120 


AH 900. 

BH 6. 

HORIZONTAL WAVELENGTH, KM 

HLt= AH + BHtCH 

DPH1I=(90. -ABS(CLAT )/0.017453293)**2 

DHGT = 0.22 + 0.00258*(SQRT(ABS(CH)**3)) 
IF (DHGT.GT.5.) DHGT = 5. 

VDS = (11.0 ~- 2.102E-4*DPHI ) *OHGT 


VTS = (3.0 + 5.146E-4*DPHI )*DHGT 

VUS = (6.2 - 3.615E-4*DPHI )*DHGT 

VDL = (20.7 - 1.3466£-3*DPHI ) *DHGT 
VTL = 7.3°DHGT 

VUL = (31.2 - 3.503E-3*DPHI ) *OHGT 


HLS = 20. + .0125*CH*tCH 
IF(HLS.GT.400.) HLS = 400. 
HLS = (DX/HLS)**2 

HLL = (DX/HLL)**2 
RDS=SORT(HLS+(0Z/VDOS)**2) 
IF(RDOS.LE. 50.)GO TO 10 
ROS=O. 

GO Ta 20 

ROS = CORREL(RDS) 
RTS=SQRT(HLS+(DZ/VTS)**2) 
IF(RTS.LE. 50.)GO TO 30 
RTS=O. 

GO TO 40 

RTS = CORREL(RTS) 
RVS=SQRT(HLS+(DZ/VUS )**2) 
IF(RVS.LE. 50.)GO TO 50 
RVS=O. 

GO TO 60 

RVS = CORREL(RVS) 

ROL =SQRT(HLL+(DZ/VDL)**2) 
IF(ROL.LE. 50.)GO TO 70 
RDL=0. 
GO TO 80 

RDL = CORREL(RDL) 
RTL=SORT(HLL+(DZ/VTL)**2) 
IF(RTL.LE. 50.)GO TO 90 
RTL=O. 
GO TO 100 
RTL = CORREL(RTL) 
RVL=SQRT(HLL+(O0Z/VUL )**2) 
IF(RVL.LE. 50.)GO TO 110 
RVL=O. 
GO TO 120 

RVL = CORREL(RVL) 

CONT INUE 


CALL CORLAT(AS,BS.CS.DS,ES,FS.GS,HS,AIS,AUS,AKS,SPS1,SPS2,SDS1, 
1 SDS2,S7TS1,STS2.$US1,SUS2.$V51,5V52,UDS1,UDS2,V0S1,VDS2,RDS,RITS, 
2RVS) 

CALL CORLAT(AL,BL,CL,OL,EL,FL,GL,HL,AIL,AJL,AKL,SPL1,SPL2,S0L1, 
4 SDL2,STL1,STL2,SUL1,SUL2,SVL1,SVL2,UDL1,UDL2,VDL1,VDL2, 
2RDL,RTL,RVL) 
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IF(SW1 .LE. O.)SW1=0.001 
IF(SW2 .LE. O.)SW2=0.001 
AW=RVS*tSW2/SWi 
BW=SW2'SQRT(1.-RVS*RVS) 

72 = RANDOM(L ) 

ZD PPND(2Z2,L) 

Z2 RANDOM(L } 

ZT = PPND(22,L) 
DS2=AS*tDS1+BS*ZD 
TS2=CStTS1+OS*DS2+ES*2ZT 
PS2=DS24+TS2 

Z2 = RANDOM(L) 

ZO = PPND(Z2,L) 

Z2 = RANDOM(L) 

ZT = PPND(Z2,L) 
US2=FS*US1+GS*DS2+HS*ZD 
VS2=AIS*VS14AUS*DS24+AKS*ZT 

22 = RANDOM(L) 


ti 


| 


27D = PPND(Z22,L) 
Z2 = RANDOM(L) 
ZT = PPND(22,L) 


DOL2=AL*DL1+BL*Z2D 
TL2=CL*TLI+DL *DL2+EL*ZT 
PL2=DL2+TL2 

22 = RANDOM(L ) 

ZD = PPND(22.L) 

22 = RANDOM(L) 

ZT = PPND(Z2,L) 
UL2=FL*UL 14+GL*OL2+HL*ZD 
VL2=AIL*VL1+AJUL*DL2+AKL *2T 
Z2=RANDOM(L ) 
ZO=PPND(272,L) 
W2=AW*WTt+BWrzD 
P2=PS2+PL2 
D2=DS2+0L2 
T2=TS2+TL2 
IF(P2.LT.-0.9) P2= -0.9 
IF(D2.LT.-0.9) D2= -0.9 
IF(T2.LT.-0.9) 12= -0.9 
U2=US24+UL2 
V2=VS24+VL2 
UDL 1=UDL2 
UNS 1=UDS2 
VDL 1=VOL2 
VOS1=V052 
RETURN 
END 


SUBROUTINE PHASE(D1,X1,0D2,X2,0,Xx) 


PER = 870. 

IF (X2=x 1) 20,.10;20 
D= 01 

RETURN 

DA = Dt 

DB = D2 
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PER2 = PER/2. 
IF(ABS(DB-DA).LE.PER2)GO TO 30 
IF (DA.LT.PER2) DA = DA + PER 
IF (OB8.LT.PER2) DB = DB + PER 
DA = DA + (DB - DA)*(X - X1)/(X2 - X1) 
IF (DA.GT.PER) DA = DA - PER 
IF(DA.LT.O. )DA=DA+PER 
D = DA 
RETURN 
END 
DOUBLE PRECISION FUNCTION PPND(P, IFAULT) 


ALGORITHM AS 114 APPL. STATIST. (1977) VOL. 26, P. 118 


PRODUCES NORMAL DEVIATE CORRESPONDING TO LOWER TATL AREA OF P. 


RETURNS IFAULT = 1 IN INPUT P >= 1 OR <= O, IFAULT = 0 
OTHERWISE. IF TFAULT = 1, PPND VALUE IS SET TO O. 

SINGLE PRECISION VERSION WITH ERROR EPSILON = 2 ** eee it 

FOR DOUBLE PRECISION VERSION, CHANGE REAL TO DOUBLE PRECISION 
IN THE FUNCTION STATEMENT AND THE DECLARATION OF VARIABLES; 
CHANGE EO TO OO IN THE DATA STATEMENTS AND CHANGE ABS, ALOG 


AND SQRT TO DABS, DLOG AND DSQRT IN THE ASSIGNMENT STATEMENTS. 


THE HASH SUMS ARE THE SUMS OF THE MODULI OF THE COEFFICIENTS. 
THEY HAVE NO INHERENT MEANINGS, BUT ARE INCLUDED FOR USE IN 
CHECKING TRANSPOSITIONS. 


DOUBLE PRECISION ZERO, SPLIT, HALF, ONE, AO, At, A2, AS, 
Bi, B2, B83, B4, CO, C1, C2, C3, DI, D2, P, Q, R 


DATA ZERO, HALF, ONE, SPLIT /O.0D0, 0.500, 1.0D0, 0.42D0/ 


DATA AO / 2.50662823884D0/, 
Ai / -18.61500062529D0/, 
A2 / 41.39119773534D0/. 
A3 / -25 .44106049637D0/, 
Bi / -8.47351093090D0/ , 
B2 / 23.0833674374300/, 
B3 / ~24.06224101826D0/, 
B4 / 3. 13082909833D0/ 


HASH SUM FOR A & B = 143.70383558076 


DATA CO / -9.78718931138D0/, 
ci / -2.29796479134D0/ , 
C2. ./ 4.85014127135D0/. 
c3 / 2.32121276858D0/ , 
D1 / 3.54388924762D0/, 
D2 / 1.6370678 1897D0/ 


HASH SUM FOR C & D = 17.43746520924 


IFAULT = O 
Q = P - HALF 
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IF(DABS(Q) .GT. SPLIT) GOTO 1 


R=Q+*9Q 
PPND = Q * (((A3 * R + A2) * R + At) * R + AO) / 
& ((((B4 * R + B39) * R + B2) * R + Bt) * R + ONE) 
RETURN 
R= Pp 


IF (Q .GT.ZERO) R = ONE - P 

IF (R .LE.ZERO) GOTO 2 

R =DSORT(-DLOG(R)) 

PPND = (((C3 * R + C2) + R + C1) * R + CO) / 
& ((02 * R + O01) * R + ONE) 

IF (Q .LT. ZERO) PPND = -PPND 

RETURN 

IFAULT = 14 

PPND = ZERO 

RETURN 

. END 

SUBROUTINE QBOGEN 
.COMPUTES QBO VALUES PQ,00.1TQ,U0,VO AT HEIGHT H, LATITUDE PHI 

ON JULIAN DAY XMJD FROM ARRAYS OF AMPLITUDES PAQ,DAO, TAQ, 
UAQ,VAQ AND PHASES PDOQ,DDO,TNO,UDO.VvOa. 

COMMON /IOTEMP/ ITOTEM1, IOTEM2,1US,DDD,XMJD,PHI1,PHI, 
NSAME,RP1, RD1, RT1, SPt, SD1, ST1, RU, RVI. SUT. 
$ MN, IDA, IYR, H1, PHIIR, THETIR.G,RI,H.PHIR, THETR.F10.F10B. AP, 
& IHR,MIN,NMORE ,DX,HL,VL,DZ,DUMMY2(24) 

COMMON /PDOTCOM/IU4,MONTH,IOPR, 

PSP(15,19,18),0SP(15,19,18),17SP(15,19,18),.USP(15,19, 18), 
VSP(15,19,18), 
PG(21,19),0DG(21,19),17TG(21,19),UG(21,19), 
PAQ(17,5),DAQ(17,.5),TAQ(17,5),UAQ(17.5),VAQ(17.5). 
PDQ(17,5),000(17,.5),TDO(17,5),UD0(17,5),VDQ( 17.5), 
PR(29,19),0R(29,19),TR(29,19),UR(29,19),VR(29, 19) 

. ,PQ.0Q,T0,U0.VO 
& ,PA,DA,TA,UA,VA, IOPQ, DUMMY (2250) 

IF (XMJD.GT.O.AND.IOPQ.EQ.1) GO TO 10 

SETS QBO VALUES TO ZERO FOR ANNUAL MEAN 

PQ =O. 

DOQ=O. 

TQ=O. 

UQ=O. 

VO=O0. 

RETURN 

LOWER HEIGHT INDEX 

IH = INT((H-5.)/5. ) 

IF (IH.LT.1) IH=1 

UPPER HEIGHT INDEX 

Ip = JH + 1 

IF (IP.GT.17) IP = 17 

PHA = ABS(PHI ) 

LOWER LATITUDE INDEX 


JL = INT(( PHA + 10.)/20. } 
UPPER LATITUDE INDEX 
JP = JL + 1 
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IF (JL.LE.O) JL=1 

IF (JP.G1T.5) UJUP=5 

JULIAN DAY FOR JAN O, 1966 

XMJDO = 2439126 

TIME RELATIVE TO JAN O, 1966 

TMJD = XMJD-XMUDO 

2*PI/PERIOD,PERIOD = 870 DAYS 

PER = 870. 

TP = 6.2831853/PER 

LOWER HEIGHT 

HI = 5. 4 §.*TH 

LOWER LATITUDE 

PHIU = 20.*JL - 10. 

UPPER LATITUDE 

PHIP = 20.*JP-10. 

INTERPOLATES QBO P,D0,T AMPLITUDE ON LATITUDE AT LOWER HEIGHT 
CALL INTERZ(PAQ( IH, JL), OAQ( IH, JL), TAQ( IH, JL), PHIU,PAQ( IH, UP), 


{DAQ(IH, JP), TAQ( IH, JP),PHIP,PA1,DA1,TA1,PHA) 
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UPPER HEIGHT 

HP = 5.+5.*IP 

INTERPOLATES OBO P,D,T AMPLITUDE ON LATITUDE AT UPPER HEIGHT 
CALL INTERZ(PAQ(IP,JUL),DAQ(IP, JUL). TAQ(C IP, JL). PHIJU,PAQ(IP,UJP), 


2DAQ(IP,JP),TAQCIP,JP),PHIP,PA2,DA2,TA2,PHA) 
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INTERPOLATES QBO P,0,T AMPLITUDE ON HEIGHT AT LATITUDE PHI 
CALL INTERZ(PA1,DA1,TA1,HI,PA2,DA2,TA2.HP,PA,DA,TA,H) 
INTERPOLATES QBO P,0D,T,U,V PHASE ON LATITUDE AND HEIGHT 

CALL PHASE(PDQ(IH,JUL),PHIJU,PDQ(IH,JP),PHIP,PD1,PHA) 

CALL PHASE(DDO( IH, JL). PHIU,ODQ( IH, JP),PHIP.DD1, PHA) 

CALL PHASE(TDQ(IH,JL),PHIU,TOQ(C IH. JP). PHIP,TOt, PHA) 

CALL PHASE(POQ(IP,JL),PHIU.POQ(IP,UP).PHIP,.PD2,PHA) 

CALL PHASE(DDQ(IP,JL),.PHIU,ODQ(IP,JP),PHIP,DD2,PHA) 

CALL PHASE(TDQCIP,JUL),PHIU, TDQ(IP,JP),PHIP,1TD2,PHA) 

CALL PHASE(PD1,HI,PD02,HP,PD,H) 

CALL PHASE(DD1.HI,DD2,HP,DD,H) 

CALL PHASE(TD1,HI,TD2,HP,TD,H) 

CALL PHASE(UDQ(IH,JL),PHIU,UDO(IH,JP),PHIP,UD1, PHA) 

CALL PHASE(VDO(IH,JL),.PHIU, VDQ( IH, JP),PHIP,VD1,PHA) 

CALL PHASE(UDQCIP,JL),FPHIU,UDQCIP,JP),PHIP,UN2,PHA) 

CALL PHASE(VDQ(IP,JL),PHIU,VDOCIP,JP),PHIP,VD2,PHA) 

CALL PHASE(UD1,HI,UD2,HF,UD,H) 

CALL PHASE(VD1,HI,VD2,HP,VD,H) 

INTERPOLATES QBO WIND AMPLITUDE ON LATITUDE AT LOWER HEIGHT 
CALL INTERWC(UAQC IH, JL), VAQ( IH, JL), PHIU,UAQ(IH, JP),VAQ(IH,JUP), 


5PHIP,UA1,VA1,PHA) 
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INTERPOLATES QBO WIND AMPLITUDES ON LATISUDE AT UPPER HEIGHT 
CALL INTERW(UAQ(IP,JL),VAQ(IP,JL),PHIU,UAQCIP,JP),VAQCIP,UP), 


GPHIP,UA2,VA2,PHA) 
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INTERPOLATES QBO WIND AMPLITUDES ON HEIGHT AT LATITUDE PHI 
CALL INTERW(UA1,VA1,HI,UA2,VA2,HP,UA,VA,H) 

EVALUATES OBO VALUES FROM INTERPOLATED AMPLITUDES AND PHASES 
PQ=PA*COS(TP*(TMUD-PD) ) 

DQ=DA*COS(TP*( 1IMJD-DD) ) 

TQ=TA*COS(TP*(TMJD-TD) ) 

UQ=UA*COS(TP*(TMJD-UD) ) 
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VQ=VA*COS(TP*(€TMJD-VD) ) 


RETURN 
END 
DOUBLE PRECISION FUNCTION RANDOM(L ) 
C 
C ALGORITHM AS 183 APPL. STATIST. (1982) VOL. 31, P.188 
C 
C RETURNS A PSEUDO-RANDOM NUMBER RECTANGULARLY DISTRIBUTED 
C BETWEEN O AND 1. 
c 
c IX, IY AND IZ SHOULD BE SET TO INTEGER VALUES BETWEEN 
C 1 AND 30,000 BEFORE FIRST ENTRY. 
c 
C INTEGER ARITHMETIC UP TO 30323 IS REQUIRED. 
C 
C RETURNS t = O UNLESS RANDOM = O OR RANDOM = 1, IN WHICH 
c CASE L = 1 
C 
DOUBLE PRECISION ONE,ZERO 
COMMON /RAND/ IX, IY, 12 
DATA ONE ,Z7ERO/1.000,0.000/ 
IX = 171 * MOD(IX, 177) - 2 * (1X / 177) 
IY = 172 * MOO(IY, 76) - 35 * (IY / 176) 
IZ = 170 * MOD(IZ, 178) - 63 * (12 / 178) 
C 
IF (IX .LT. O) IX = IX + 30269 
IF (IY .LT. O) TY = LY + 30307 
IF (IZ .LT. O) 12 = 12 + 30323 
c 
C IF INTEGER ARITHMETIC UP TO 5,212,632 IS AVAILABLE, 
C THE PRECEDING 6 STATEMENTS MAY BE REPLACED BY 
C 
Cc IX = MOD(171 * IX, 30269) 
C IY = MOD(172 * IY, 30307) 
C IZ = MOO(170 * 12, 30323) 
C 
C ON SOME MACHINES, THIS MAY SLIGHTLY INCREASE THE SPEED. 
Cc THE RESULTS SHOULD BE IDENTICAL. 
C 
RANDOM =DMOD(DFLOAT( IX) /30269.O0D0+DFLOAT(IY) /30307.0D0+ 
& DFLOAT(IZ) /30323.OD0, ONE ) 
ix = 6 
IF (RANDOM .LE. ZERO .OR. RANDOM .GE. ONE)L = 1 
RETURN 
ENO 
SUBROUTINE RIG 
COMMON/JTOTEMP/IOTEM1, IOTEM2,1US,00.XMJD,PHI41,PHI, 
NSAME ,RP1, RDI, RT1, SP1, SD1, ST1, RUT, RV1, SU1, 
¢ MN, IDA, IYR, H1, PHIIR, THETIR,G,RI,H,PHIR, THETR,F10,F 108, AP, 
$ IHR .MIN,NMORE ,DX,HL,VL,DZ,8,EPS, IOPP,LOOK,IET,GLAT, 
IRP1S,RD1S,RT1S,RUIS,RV1S,5P15,SD1S,S5T1S,SU1S,SV1S, 
2UDS1,VDSt,UDL1,VDL1,UDS2.VD0S2,UDL2,VDI_2 
Cy. 5 GRAVITY G AT H, LATITUDE PHIR (RADIANS) 
C..... RADIUS RI FROM CENTER OF EARTH TO HEIGHT H 
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B = POLAR EARTH RADIUS, EPS = ECCENTRICITY 
CPHI2 = COS(PHIR) ** 2 
EARTH RADIUS 
RI = B / SQRT(1. - EPS * CPHI2) 
C2PHI = COS(2*PHIR) 
C2PHI = 2. * CPHI2 - 1. 
C4PHI = COS(4*PHIR) 
C4PHI = 8. * CPHI2 * (CPHI2 - 1.) + 1. 
G AT SURFACE 
G = 9.80616 * (1. - 0.0026373 * C2PIII + 0.0000059 * C2PHI 
EFFECTIVE RADIUS 


RE = 2. * G / (3.085462E-3 + C2PHI * 2.27F-G - C4PHI * 2.E-9) 


G AT HEIGHT H 
G= G/ (1. + (H / RE)) ** 2 
RADIUS AT HEIGHT H 
RI = RI + H 
END 
SUBROUTINE RTERP(H,PHI,PR,DR,.TR,P,D,T) 
COMPUTES RANDOM PERTURBATION STANDARD DEVIATIONS P,D.T AT 
HEIGHT H (KM), LATITUDE PHI(DEGREES) FROM SIGMA ARRAYS 
PR,DR,AND TR 
DIMENSION PR(29,19),0R(29,19),TR( 29, 19) 
I = LOWER HEIGHT INDEX 
IF(H .LT. 120. )I=INTC(CHtS. )/5.) 
IF (tl .GE. 120. )1=25+INT((H-1420. )/20. ) 
IF(T.UT.1)2 = 1 
IF(I.GT.29)I = 29 
IP = I#14 
IF(IP .GT. 29)IP=29 
LOWER LATITUDE INDEX 
J=INTC (PHT + 100.)/10.) 
JP = Jif 
IF(JP .GT. 19)JP=19 
IF(1l1 .GT. 25)GOTO 10 
LOWER HEIGHT FOR PR,TR,DR ARRAYS 
ZA25 6" 1S. 
GO TO 20 
Z71=120.+20.*( I-25) 
IF(IP .GT. 25)GOTO 30 
UPPER HEIGHT FOR PR,DR,TR ARRAYS 
Z2=5.+*1P-5. 
GO TO AO 
Z2=120.+20.*CIP-25) 
PHI1=-100.410.tJ 
PHI2=-100.+10. *UP 
INTERPOLATE ON LATITUDE AT LOWER HEIGHT 


CALL INTERZ(PROI,J) ORCI, J), TROL, J). PHI, PROT, JP),DR(™ JP), 


TREY PY PHTID PI M4, Tt. PHT) 
INTERPOLATE ON LATITUDE AT UPPER HEIGHT 


CALL INTERZCPRCIP a) ,OREIPL J). TROEIPL J), PHIT,PROEIP, JP) ,DORCIP. JP), 


TR( IP, JP), PHI2,P2,02,T2. PHT) 
INTERPOLATION ON HETGHT USING LATITUDE INTERPOLATED VALUES 
CALL INTFRZ(P1,01,71.21,P2.02,7T2,22,P,0,1T.H) 
RETURN 
END 
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SUBROUTINE RTRAN RTRA 
COMMON/IOTEMP/IOTEM1, IOTEM2,1IUS,DOUMMY(6Q) RTRA 
CHARACTER*2 I1 RTRA 
COMMON/COTRAN/NDATA(19).12,13,15,14(10) RTRA 
COMMON/COTRAN1/I1 RTRA 

C.... ENTRY POINT TO READ STATIONARY PERTURBATION DATA RTRA 
CHARACTER IWHERE*10 RTRA 
IWHERE=‘RTRAN. RTRA 
READ(IUS,10,END=3) I11,12,I13,NDATA RTRA 

10 FORMAT(A2,2115) RTRA 
Is = 0 RTRA 
RETURN RTRA 

ENTRY RTRAN‘1 RTRA 

C.... ENTRY POINT TO READ ZONAL MEAN DATA & RANDOM PERTURBATION DATA IN RTRA 
Cc SETUP RTRA 
IWHERE=‘RTRAN1. RTRA 
READ(IUS,200,END=3) I1,.12,I13,NDATA,1I5 RTRA 

200 FORMAT(A2.14,15,1916,14) RTRA 
RETURN RTRA 

ENTRY RTRAN2 RTRA 

C.... ENTRY POINT TO READ QBO PARAMETERS AND LARGE/SMALL SCALE FRACTIONSRTRA 
C.... AND CORRELATIONS IN SETUP RTRA 
IWHERE=‘RTRAN2. RTRA 
READ(IUS,300,END=3) I1.NDATA RTRA 

300 FORMAT(A2, 1917) RTRA 
I2=NDATA(1) RTRA 
I3=NDATA(2) RTRA 

15=0 RTRA 

DO 2 I = 1,10 RTRA 

2 I4(1)} = NDATA(1+I1) RTRA 
RETURN RTRA 

3 WRITE(6,250) IUS, IWHERE RTRA 
250 FORMAT(’1 PREMATURE END-OF-FILE FOUND ON UNIT “’,12/ RTRA 
1’0 CALLED FROM SUBROUTINE ‘’,A10) RTRA 

STOP ’ RTRAN ERROR.‘ RTRA 

END RTRA 
SUBROUTINE SCIMOD(NPOP ) SCIM 
Ce COMPUTES VALUES P,D,T,U.V AND SHEAR DUH,DVH FROM INPUT AND SCIM 
Cc ARRAYS IN COMMON POTCOM. INPUT TO SCIMOD IS% SCIM 
C G = GRAVITY AT POSITION RI = RADIUS AT HEIGHT H SCIM 
Cc PHIR = LATITUDE (RADIANS) THETR = LONGITUDE (RADIANS) SCIM 
Cc FIO = F10.7 SOLAR FLUX FIOB = MEAN F10.7 FLUX SCIM 
Cc AP = SOLAR-GEOMAGNETIC A SUB P INDEX SCIM 
Cc MN/IDA/IYR = DATA (IYR = FULL YEAR-1900) SCIM 
C IHRYMIN = TIME H1 = PREVIOUS HEIGHT SCIM 
Cc PHI1R = PREVIOUS LATITUDE THETIR = PREVIOUS LONGITUDE SCIM 
c RP1,RO1,RT1 = PREVIOUS RANDOM PERTURBATIONS SCIM 
C SP1,5D1,.ST1 = PREVIOUS RANDOM STANDARD DEVIATIONS (SIGMAS) SCIM 
C RU1,RVi = PREVIOUS RANDOM WINDS SCIM 
C SU1,SV1 = PREVIOUS RANDOM WIND SIGMAS SCIM 
COMMON/IPRTP/ IPRI,NLIMIT SCIM 
COMMON/IOTEMP/IOTEM1, IOTEM2,IUS,0D,XMUD,PHI1,PHI, SCIM 
-NSAME ,RPIL ,RDIL.RTIL,SPIL, SDIL,STIL,RUIL,RVIL,SUIL,SVIL. SCIM 

¢ MN. IDA, TYR, Ht, PHIIR, THETIR,G,RI,H,PHIR, THETR,F10.F 108, AP. SCIM 
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IHR,MIN,NMORE ,DX,HL,VL,0Z,B,EPS,IOPP,LOOK,IET,FLAT, 
{RP1S,RDIS,RT1IS,RUIS,RV1IS,SP1S,5D15.5T15,S5U1S,5V15S, 
2UDS1,V0S1,UDLt,VDL1,U0S2,VDS2,UDL2,VDL2 

COMMON /PDTCOM/1U4,MONTH,IOPR, 

. PSP(15,19,18),0SP(15, 19.18), TSP(15,19,.18),USP(15,19,18), 


PG(21,19),0G(21,19),TG(21, 19) ,UG(21, 19), 
PAQ(17.5),.0AQ(17.5),TAQ( 17.5), UAQ( 17.5), VAQ( 17,5), 
PDQ(17,5),000(17,5),.10Q(17,5),UD0(17,5),VD0( 17,5), 
PR(29,19),0R(29,.19),TR(29,19),UR(29, 19), VR(29,19) 
. .PQ,00,TO,UQ, VO,POA,00A,TOA,UA,VA,IOPO, 
1PLP(25,10),OLP(25,10),TLP(25, 10), 
2ULP(25,10),VLP(25,10),UDL(25, 10), 
3VDL(25,10),UDS(25, 10), VDS(25, 10) 
COMMON /C4/ GLAT(16),GLON(16),NG,P40(16,26),D4D(16,26),140( 16,26), 
SP4(16,26),504(16,26).S5T4( 16,26), THET1, THET , DUMMY 
COMMON/COMPER/SPH, SOH, STH, PRH, DRH, TRH, URH, VRH, SUH, SVH,CP, 
4PRHS ,ORHS, TRHS,URHS, VRHS,PRHL ,DRHL. TRHL,URHL , VRHL, 
2SPHS , SOHS. STHS, SUHS, SVHS, SPHL , SOHL , STHL , SUHL , SVHL 
COMMON/WINCOM/DH.FCORY ,DX5,DY5,0PX,DPY ,UGH, VGH, 
TH ,OTX ,DTY ,DUH, DVH,PH,UPRE , VPRE , DUPRE , DVPRE 
COMMON/CHK/PCK(4,4,3),0C0K(4,4,3),NO(2) 
COMMON /CHIC/LA(4,4),NB(2), IWSYM,USH, VSH,DUSH,DVSH 
COMMON /VERT/RW1,SW1,WRH, SWH, WR(29) 
FACTOR FOR RADIANS TO DEGREES 


FAC = 


IWSYM=ICHAR( ’ 


VSP( 15,19, 18), 


57.2957795 


#) 


IF(NPOP.NE.O) GO TO & 
UPRE=0O. 
VPRE=O. 


DUPRE= 
DVPRE= 


PQ=O. 
DO=O0. 
TQ=0. 
PRH=0, 


DRH=0. 
TRH=0. 
URHI=O. 


VRH=O. 
WRH = 
UQ=0. 
VQ=0. 
PQA=0. 


DQA=0. 


TQA=0. 
UA=0. 
VA=0. 
PSH=O0. 
DSH=O. 
TSH=0. 
SPU = 
SPV = 


QO. 
0. 


0. 
0. 
MONTH=MN 
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PRESENT LATITUDE, DEG 

PHI = PHIR*FAC 

PRESENT LONGITUDE, DEG 

THET = THETR*FAC 

PREVIOUS LATITUDE, DEG 

PH11 = PHIIR*FAC 

PREVIOUS LONGITUDE, DEG 

THET1 = THETIR*FAC 

FCORY = NORTH COMPONENT CORIOLIS FACTOR TIMES DISTANCE FOR 
5S DEGREES OF LATITUDE 

DYS5 = 5000.*RI/FAC 

DX5 = DYS*COS(PHIR) 

FCORY = DYS*SIN(PHIR)/( 120. *FAC) 

IN JACCHIA OR MIXED ZONAL MEAN-JACCHIA HEJGHT RANGE 

IF(H.GT.90.0) GO TO 10 

IN 4-D DATA HEIGHT RANGE 

IF (H.LE.25.0) GO TO 500 

IN ZONAL MEAN OR MIXED ZONAL MEAN 40 HEIGHT RANGE 

GO TO 200 

IN MIXED JACCHIA-ZONAL MEAN RANGE, NEED TO FAIR DATA 

IF (H.LT.120.) GO TO 20 

FOLLOWING IS THE PURE JACCHIA HEIGHT RANGE SECTION 

JACCHIA VALUES AT CURRENT POSITION 

CALL JACCH(H,PHIR, THET,PH,DH, TH) 

PHIN = PHIR + 5S. / FAC 

THETE = THET - 5. 

JACCHIA VALUES AT CURRENT POSITION+5S DEGREES LAT, FOR DP/DY AND 
DT/DY 

CALL JACCH(H,PHIN, THET,PHN,DHN, THN) 

JACCHIA VALUES AT CURRENT POSITION-5 DEGREES LON, FOR DP/DX AND 
DT /DX 

CALL JACCH(H,PHIR, THETE ,PHE,.DHE. THE ) 

DP/DY FOR GEOSTROPHIC WIND 

DP Y=PHN-PH 

OP/DX FOR GEOSTROPHIC WIND 

DPX =PHE-PH 

DT/DX FOR THERMAL WIND SHEAR 

OTX = THE - TH 

DT/DY FOR THERMAL WIND SHEAR 

DTY = THN - TH 

CHANGE NOTATION FOR OUTPUT 


PGH=PH 

DGH=DH 

1 GH=TH 

CALL WIND 

UH = UGH 

VH = VGH 

HB = H+ 5. 

CP = 7.*PH/(2.*DH*TH) 


CALL JACCH(HB,PHIR, THET,PB,0B,TB) 
OTZ = (TB - TH)/SOOO. 
VERTICAL MEAN WIND 


WGH = -CP*(UH*DTX/DX5 + VH*DTY/DY5)/(G + CP*DTZ + UH*DUH+VH*DVH) 


GO TO RANDOM PERTURBATIONS SECTION 
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GO TO 800 


FOLLOWING IS THE MIXED JACCHIA-ZONAL MEAN HEIGHT RANGE SECTION 


LOWER HEIGHT INDEX 
IHA = S*CINT(H)/5) 
UPPER HEIGHT INDEX 
IHB = IHA + § 
LOWER HEIGHT FOR INTERPOLATION 


HA = ITHA*1. 
UPPER HEIGHT FOR INTERPOLATION 
HB = IHB*1. 


JACCHIA VALUES AT LOWER HEIGHT, CURRENT LAT-LON 
CALL JACCH(HA,PHIR, THET,PJA,OJUA.TJA) 
PHIN = PHIR + 5. / FAC 
THETE = THET - 5. 
JACCHIA VALUES AT LOWER HEIGHT, CURRENT LON-LAT+5 DEGREES 
LAT, FOR DP/DY AND OT/DY 
CALL JACCH(HA,PHIN, THET, PUN, DUN, TUN) 
JACCHIA VALUES AT LOWER HEIGHT, CURRENT LAT-LON-5 DEGREES 
LON, FOR DP/DX, AND DT/0x 
CALL JACCH(HA,PHIR, THETE,PUJE,DUJE, TJE) 
JACCHIA DP/DY AT LOWER HEIGHT 
DPXJA=PJE-PJA 
JACCHIA DP/DY AT LOWER HEIGHT 
DPYJA=PUN-PJA | 
JACCHIA DT/DX AT LOWER HEIGHT 
DTXJA = TJE - TJA 
JACCHIA DT/DY AT LOWER HEIGHT 
OTYJA = TJN - TJA 
JACCHIA VALUES AT UPPER HEIGHT, CURRENT LAT-LON 
CALL JACCH(HB,PHIR, THET,PJB,DUB, TUB) 
PHIN = PHIR + 5. / FAC 
THETE=THET-5. 
JACCHIA VALUES AT UPPER HEIGHT, CURRENT LON-LAT+5 DEGREES 
LAT, FOR DP/DY AND Di/DbDyY 
CALL JACCH(HB,PHIN, THET,PUJUN,DUN, TUN) 
JACCHIA VALUES AT UPPER HEIGHT, CURRENT LAT-LON-5 DEGREES 
LON, FOR DP/DX AND DT/DXx 
CALL JACCH(HB,PHIR, THETE ,PJUE,DJE, TUE) 
JACCHIA DP/DX FOR GEOSTROPHIC WINDS 
DPXJB = PJE - PJB 
JACCHIA DP/DY FOR GEOSTOPHIC WINDS 
DPYJB = PUN ~- PUB 
VYJACCHIA DT/DX FOR THERMAL WIND SHEAR 
DTIXJB = TJE - TUB 
JACCHIA DT/DY FOR THERMAL WIND SHEAR 
DTYJB = TUN - TuJB 
ZONAL MEAN AT LOWER HEIGHT, TO BE FAIRED WITH JACCHIA 


CALL GTERP(IHA,PHI,PGA,DGA,TGA,PG,DG,TG,DPYGA,DTYGA,UGA,UG) 


ZONAL MEAN AT UPPER HEIGHT, TO BE FAIRED WITH JACCHIA 


CALL GTERP(IHB,PHI,PGB,DGB,TGB,PG,DG, TG,DPYGB,DTYGB,UGB,UG) 


FATRED RESULTS AT LOWER HEIGHT 
IHSB = 90 
CALL PDTUV(PSP,DSP,TSP,USP,VSP,PHI, THET, IHSB,PSH,OSH, TSH, 


& DPXSB,DPYSB, DTXSB,DTYSB,SPU,SPV) 
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PGA = PGA*t(1. + PSH) SCIM 

DGA = DGAt(1t. + DSH) SCIM 

TGA = TGA*(1. + TSH) SCIM 

PGB = PGB*(1. + PSH) SCIM 

DGB = DGBt(1. + DSH) SCIM 

TGB = TGB*(1. + TSH) SCIM 
UGA=UGA+S5PU SCIM 
VGA=SPV SCIM 

UGB =UGB+SPU SCIM 
VGB=SPV SCIM 
DTXGA = DTXSB * TGA SCIM 
DTXGB = DTXSB * TGB SCIM 
DTYGA = TGA*DTYSB + OTYGA*(1. + TSH + OTYSB) SCIM 
DTYGB = TGB*DTYSB + OTYGB*(1. + TSH + DTYSB) SCIM 
OPXGA = DPXSB * PGA SCIM 
DPXGB = DPXSB * PGB SCIM 
DPYGA = PGA*OPYSB + DPYGA*(1. + PSH + DPYSB) SCIM 
DPYGB = PGB*DPYSB + DPYGB*(1. + PSH + DPYSB) SCIM 

CALL FAIR(PGA,.DGA,TGA.PJA,OJUA,TJA,IHA,.P1,.D1,1T1,DPXGA,DPYGA, SCIM 

& OPXJA,DPYJUA,DPXA,DPYA,DTXGA,DTYGA,DTXJUJA,DTYJA,DTXA,DTYA) SCIM 

Asta sons FAITRED RESULTS AT UPPER HEIGHT SCIM 
CALL FAIR(PGB.DGB, TGB.PJB.DUB, TUB, IHB,P2,D2.T2,0PXGB,DPYGB, SCIM 

&® DPXJB,DPYUB,DPXB,DPYB,DTXGB ,OTYGB ,OTXJUB,OTYJ8 ,DTXB,DTYB) SCIM 
aes HEIGHT INTERPOLATION ON FAIRED P,O,T SCIM 
CALL INTER2(P1,.01,1T1,HA,P2,02,172,HB,PH,DH,.TH,H) SCIM 

Ohad HEIGHT INTERPOLATION ON FAIRED DP/DX,DP/DY SCIM 
CALL INTERW(DPXA,DPYA.HA,OPXB,OPYB,HB,DPX,DPY,H) SCIM 

oh ey HEIGHT INTERPOLATION ON FAIRED DT/DX ,DT/NY SCIM 
CALL INTERW(OTXA,OTYA,HA,DTXB ,DTYB,HB,OTX,OTY,H) SCIM 
HEIGHT INTERPOLATION OF WIND SCIM 
CALL INTERW(UGA, VGA,HA,UGB, VGB,HB.USH, VSH,H) SCIM 
DUSH=(UGB-UGA)/5000. SCIM 
DVSH=(VGB-VGA)/5000. SCIM 
CHANGE OF VARIABLES FOR OUTPUT SCIM 
PGH=PH SCIM 
DGH=0H SCIM 
TGH=TH SCIM 
CALL WIND SCIM 
UH=UGH SCIM 
VH=VGH SCIM 

CP = 7.*PH/(2.*DH*TH) SCIM 

DTZ = (T2 - Tt)/5000. SCIM 
Stes VERTICAL MEAN WINO SCIM 
WGH = -CP*(UH*DTX/DX5 + VH*DTY/DYS)/(G + CP*DTZ + UH*DUH + VH*DVH)SCIM 

GO TO RANDOM PERTURBATIONS SECTION SCIM 

GO TO 800 SCIM 

THE FOLLOWING SECTION IS FOR ZONAL MEAN OR MIXED ZONAL MEAN 4D SCIM 
HEIGHTS SCIM 
UPPER HEIGHT INDEX SCIM 

200 IHGB = 5*(INT(H)/5) + § SCIM 
UPPER HEIGHT SCIM 

HGB = IHGB*1. SCIM 
ZONAL MEAN AT UPPER HEIGHT SCIM 

CALL GTERP(IHGB.PHI .PGB,OGB.TGB,PG,.DG,TG, DPYGB,DTYGB ,UGB,UG) SCIM 
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THSB = S*CINT(H)/5) + & 

IF (ITHSB .GT. SO)IHSB = 30 

UPPER STATIONARY PERTURBATION HEIGHT 

HSB = IHSBT4. 

STATIONARY PERTURBATIONS AT UPPER HEIGHY 

CALL POTUV(PSP,OSP,TSP,USP, VSP, PHI. THET, IHS8,PSB,DSB,TSB, 


$ DPXSB,DOPYSB, DOTXSB,DTYSB,USB, VSB) 
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LOWER HEIGHT INDEX 

THGA = IHGB - 5 

LOWER HETGHT INDEX 

HGA = THGA*1. 

ZONAL MEAN AT LOWER HEIGHT 

CALL GTERP(TIHGA,PHI,PGA,DGA,TGA,PG,0G, TG, OPYGA,.DTYGA,UGA,UG) 
THSA=IHSB - 5 

LOWER STATIONARY PERTURBATION HEIGHT 

HSA = THSA*1. 

STATIONARY PERTURBATIONS AT LOWER HEIGHT 

CALL PDTUV(PSP,DSP,TSP,USP,VSP,PHI.THET, IHSA,PSA,DSA,TSA, 


$ DPXSA,OPYSA, OTXSA,DTYSA,USA,VSA) 
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CALL INTERW(UGA,O.,HGA,UGB,O. ,HGB,UGH, VGH,H) 

CALL INTERW(USA, VSA,HSA,USB, VSB,HSB, SPU, SPV,H) 

USH=UGH+4 SPU 

VSH=SPV 

DUSH=( (UGB-UGA ) / (HGB-HGA )+(USB-USA)/(HSB-HSA) )*.001 
DVSH=.001*(VSB-VSA)/(HSB-HSA) 

FOR MIXED ZONAL MEAN - 4D SECTION 

IF(H.LT.30.0) GO TO 300 

ZONAL MEAN VALUES HEIGHT INTERPOLATIONS 

CALL INTER2(PGA,OGA, TGA,HGA,PGB,DGB,TGB,HGB,PGH,OGH, TGH,H) 
STATIONARY PERTURBATION HEIGHT INTERPOLATION 

CALL INTERZ(PSA.0OSA,TSA,HSA,PSB,DSB, TSB,HSB.PSH, DSH, TSH,H) 
QUASI-BIENNIAL VALUES 

CALL QBOGEN 

HEIGHT INTERPOLATION OF ZONAL MEAN DP/DY AND DT/DY 

CALL INTERW(DPYGA .DTYGA,HGA,OPYGB ,DTYGB,HGB,DPYG, 


$ DTYG,H) 
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HEIGHT INTERPOLATION OF STATIONARY PERTURBATION DP/DX AND OP/DY 
CALL INTERW(DPXSA,DPYSA,HSA,OPXSB.DPYSB.HSB.OPXS,DPYS,H) 

HETGHT INTERPOLATION OF STATIONARY PERTURBATION DT/OX AND OT/DY 
CALL INTERW(DTXSA.OTYSA,HSA,OTXSB .DIYSB,.HSB,DTXS,DTYS,H) 
UNPERTURBED (MONTHLY MEAN) VALUES FOR OUTPUT 

TGH = TGH * (1. + TSH) 

PGH PGH * (1. + PSH) 

DGH DGH * (1. + DSH) 

TOTAL DT/DX 

DTX = DTXS + TGH 

TOTAL DT/DY 

DTY = TGHtDTYS + DTYGt(1. + TSH + DTYS) 

TOTAL DP/DxX 

DPX = DPXS * PGH 

TOTAL DP/DY 

DPY = PGH*DPYS + OPYG*(1. + PSH + DPYS) 

UNPERTURBED VALUES PLUS QBO PERTURBATIONS 

PH = (1. + PQ) * PGH 
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DH = DGH * (14. + DQ) 

TH = (4. + TQ) * TGH 

CALL WIND 

GEOSTROPHIC WIND PLUS QBO WIND PERTURBATIONS 

UH=UGH+U0 

VH=VGH+ VQ 

CP = 7.*PGH/(2.*OGH*TGH) 

OTZ = (TGBt(1.+TSB) - TGA*(t.+TSA))/5000. 
VERTICAL MEAN WIND 

WGH=-CP* (UGH*DTX/DX54+VGH*DTY/DYS)/(G+CP*DTZ+VGH* DUH+ VGH*DVH) 
GO TO RANDOM PERTURBATIONS SECTION 

GO TO 800 

THE FOLLOWING IS THE MIXED ZONAL MEAN-4D SECTION 

GENERATE GRID OF 4D PROFILES IF PREVIOUS HEIGHT GE 30 
IF (LOOK .EQ. 1)CALL GEN4D 

IHCK = 24 

DO 310 KND = 1,3 

IKND = IHCK + KNOD 

IF (IKND.GT.26)1IKND=26 

DO 310 INO 1,4 

00 310 JND 1,4 

PCK(IND.UND,KNO) = P4D(4*( IND-1)+UND,IKND) 

DCK(IND,JND,KNO) = D40(4*(IND-1)+JND, IKND) 

CONT INUE 


--LAT-LON INTERPOLATION OF 40D DATA AT 25 KM 


CALL INTER4( PHI ,THET,25, P4D ,D40,T4D,P4A,D4A,T4A, 


$ DPX4,0PY4,DTX4,DTY4) 


350 


ZONAL MEAN PLUS STATIONARY PERTURBATIONS 
PB = PGBt(1. + PSB) 
P.D,T 
DB = DGB+(1. + DSB) 
TR = TGB*(1. + TSB) 
OPXB = PGB*DPXSB 
DPYB = PGB*DPYSB + DPYGB*(1. + PSB + DPYSB) 
DTXB = TGB*DTXSB 
OTYB = TGB*DTYSB + DTYGB*(1. + TSB + DTYSB) 
HEIGHT INTERPOLATION BETWEEN 4D AT 25 AND ZONAL MEAN AT UPPER 
HEIGHT DP/DX AND DP/DY 
CALL INTERW(DPX4,O0PY4, 25. .0PXB,.O0PYB,HSB,DPX.DPY,H) 
HEIGHT INTERPOLATION BETWEEN 40 AT 25 ANDO ZONAL MEAN AT UPPER 
HEIGHT P,D,T . 
CALL INTER2(P4A,04A,T4A,25.,PB,D8, TB, HGB,PGH,DGH, TGH,H) 
HEIGHT INTERPOLATION BETWEEN 4D AT 25 AND ZONAL MEAN AT UPPER 
HEIGHT DT/DX AND DT/DY 
CALL INTERW(DTX4 .DTY4,25. ,DTXB,DTYB,.HSB.DTX ,DTY,H) 
IF (IOPQ.EQ.2) GO TO 350 
QUASI BIENNIAL PERTURBATIONS 
CALL QBOGEN | 
ADD OBO PERTURBATIONS TO P,O,T 
PH=PGH*(1.+PQ) 
DH=DGH*(1.+DQ) 
TH=TGH*(1.+TQ) 
CALL WIND 
ADD QBO WIND PERTURBATIONS 
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UH=UGH+UQ 
VH=VGH+VQ 
CP = 7.*PGH/(2.*DGH*TGH) 
DTZ = (TB - T4A)/(1000.*(HGB - 25.)) 
VERTICAL MEAN WIND 
WGH= -CP* (UGH*D1X/DX54+VGH*DTY/DY5 )/(G+CP *DT7+UGH*DUH+VGH*DVH) 
GO TO RANDOM PERTURBATIONS SECTION 
FORMAT(’ 4-D DATA AFTER ADUJUSTMENTS’/’ LATITUDE’/3X, 16F8.3) 
FORMAT(’ LONGITUDE’ /3X, 16F8.3) 
FORMAT(’ PRESSURE’ ) 
FORMAT(1X,I12,16F8.0) 
FORMAT(’ DENSITY‘ ) 
FORMAT(’ TEMPERATURE’ ) 
FORMAT(1X,12,16F8.5) 
FORMAT(1X,12,16F8.2) 
GO TO 800 
IF (H.GE.0.0) GO TO 510 
IF (H.LT.-0.015) GO TO 505 
IF ~15 METER LE HLT O ,H IS SET TOO 
H = O. 
GO TO 510 
NO MORE COMPUTATIONS TO BE MADE IF HEIGHT LT -5 M 
NMORE = O 
RETURN 
GENERATE GRIO OF 4D PROFILES IF PREVIOUS HEIGHT GE 30 
IF (LOOK .€Q. 1)CALL GEN4D 
LOWER HEIGHT INDEX 
IHA=INT(H) 
LOWER HEIGHT INDEX 
HA = JTHA*?. 
IWSX = IWSYM 
IHCK=IHA- 1 
DO 511 KND=1,3 
IKND = IHCK + KND 
IF (IKND.LT.1)IKND = 1 
IF (IKNO.GT.26)IKNOD = 26 
DO 511 IND=1,4 
DO 511 JND = 1,4 
PCK( IND, JND,KND )=P4D(4*( IND-1)+JNO, IKND) 
DCK(INO,JUND,KNO)=D4D(4*( IND-1)+JUNO, IKND) 
CONTINUE 
UPPER HEIGHT INDEX 
IHB = IHA + 4 
IF(IHB.LE.25) GO TO $13 
IHA=24 
HA=24. 
IHB=25 
UPPER HEIGHT 
HB = IHB*?. 
LAT-LON INTERPOLATION OF 4D VALUES AT UPPER HEIGHT 
CALL INTER4( PHI, THET,IHB, P4D ,D4D,T40,PB,DB.1B, 


$ DPX4B,DPY4B,DTX4B,0TY4B) 


IF(IHA.EQ.0.AND.P8*0B*TB.LE.O.)GO TO 520 
GO TO 540 
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520 IHB=IHB+1 
Benoit ae Loorp TO FIND LOWEST VALID HEIGHT 
HB=HB+1. 
GO TO 515 
540 IFC(IHA.GT.O)CALL INTER4( PHI, THET.IHA, P40 ,D4D,T4D, 
$ PA,DA,TA,DPX4A,DPY4A,DTX4A ,OTY4A) 
IFC IWSYM.EQ.ICHAR(’*’)) IWSX=IWSYM 
IF(IHA.EQ.0.0R. (PA*DA*TA.LE.O.AND.THA.LT. 10.AND.PB*DB*TB.GT.O. )) 
iGO TO 550 
GO TO 600 
ae 2 LAT-LON INTERPOLATION OF 40 VALUES AT LOWER HEIGHT 
550 CALL INTER4( PHI , THET,O, P4D.D4D.T4D. 
$ PA,DA,TA,DPX4A,DPY4A,DTX4A,DTY4A) 
IF(IWSYM.EQ.ICHAR(‘*’)) IWSX=IWSYM 
IF (TA-~1TB)560,570,560 
56Q IF(TA*TB.LE.0.0) GO TO 570 
TZ = (TA-TB) / ALOG(TA/TB) 
GO TO 575 
570 TZ=TA 
_., COMPUTES HEIGHT OF SURFACE 
575 HA = HB 
IF(PB*PA.LE.0.0)GO TO 576 
HA = HB + 0.28705*TZ*ALOG(PB/PA)/G 
576 IF(H.GT.HA - ©.04)GO TO 600 
PH=O0. 
DH=0. 
TH=O. 
PGH=0. 
DGH=0. 
TGH=0. 
GO TO 800 
ected HEIGHT INTERPOLATION OF P,D.T 
GOO CALL INTER2(PA.DA, TA,HA,PB,DB,1TB,HB,PGH,DGH, TGH,H) 
poe ane HEIGHT INTERPOLATION OF OP/DX AND DP/DY 
CALL INTERW(DPX4A,DPY4A,HA,OPX4B, DPY4B,HB,OPX ,OPY,H) 
nate HEIGHT INTERPOLATION OF DT/DX AND OT/DY 
CALL INTERWC(DTX4A,DTY4A,HA,DTX4B8 ,.DTY4B,HB,OTX ,DTY,H) 
CHANGE OF NOTATION FOR OUTPUT 
PH = PGH 
DH = DGH 
TH = TGH 
IF (PH*DH*TH.LE.O.) GO TO 800 
CALL WIND 
CHANGE OF NOTATION FOR OUTPUT 
UH = UGH 
VH = VGH 
CP = 7.*PGH/(2.*OGH*TGH) 
DTz = (TB - TA)/(1000.*(HB - HA)) 
VERTICAL MEAN WIND 
WGH = -CP*(UGHtDTX/DX5S + VGH*DTY/DY5)/(G+CP *0T Z2+UH*DUH+VH*DVH ) 
QBO=O IF HLT 10 
IF (H.LT.10.) GO TO 800 
IF (10PQ.EQ.2) GO TO 650 
COMPUTES QUASI BIENNIAL PERTURBATIONS 
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CALL QBOGEN 
ADDS QBO PERTURBATIONS TO P.D.T 
PH=PGH*(1.4+PQ) 
DH=DGH*(1.+DQ) 
TH=TGH*(1.4+TQ) 
ADDS QBO WIND PERTURBATIONS TO U,V 
UH=UGH+UQ 
VH=VGH+VQ 
THE FOLLOWING IS THE RANDOM PERTURBATIONS SECTION 
NO RANDOM PERTURBATIONS IF IOPR GT 1 
CONTINUE 
IF(H .GE. 30.)GOTO 512 
IFC IPRT.NE.O)GOTO 512 
WRITE(6,2000) (GLAT(I),1I=1,NG) 
WRITE(6,2001) (GLON{I),J]=1,NG) 
WRITE(6, 2007) 
0O 504 I=1,26 
IH=1-1 
WRITE(6,2004)1IH,(SP4(J,1),JU=1.NG) 
WRITE(6,2002) IH, (P4D(U,1I),J=1.NG) 
CONTINUE 
WRITE(6, 2003) 
DO 507 I = 1,26 
IH = 1-1 
WRITE(6, 2004) IH. (SD4(U.1),JU=1,NG) 
WRITE(6,2004)IH, (D4D(JU,1),U=1.NG) 
WRITE(6,2005) 
DO 506 I = 1,26 
IH = 1 - 1 
WRITE(6, 2004 )IH. (ST4(JU,T).JU 
WRITE(6, 2006) IH, (T4D(JU.1).J 
IPRT=IPRT+1 
CONTINUE 
IF(NPOP.EQ.0)GO TO 840 
IF (10PR.GT.1) GO TO 830 
INTERPOLATES RANDOM WIND MAGNITUDES [TO HEIGHT H, 
CALL INTRUV(UR,VR,H,PHI, SUH, SVH) 
CALL INTR25(PLP,DLP,H,PHI,PLPH, DI PH) 
CALL INTR25(TLP,DLP,H,PHI,TLPH, DLPH) 
CALL INTR25(ULP,VLP,H, PHI ,ULPH, VLPH) 
CALL INTR25(UDL,VDL,H,PHI,UDL2,VDL2) 
CALL INTR25(UDS,VOS,H,PHI,U0S2,VDS2) 
CALL INTRW(WR,H, SWH) 
SUHL =SQRT (ULPH* ABS (SUH) ) 
SUHS=SORT((1.-ULPH)*ABS( SUH) ) 
SVHL=SORT(VLPH*ABS(SVH) ) 
SVHS=SORT((14.-VLPH)*ABS(SVH)) 
SUH = SORT(ABS(SUH) ) 
SVH = SQRT(ABS(SVH)) 
IF(H.GE.25.)GOTO 805 
IF H LE 20 USE 4D DATA RANDOM P.D.T SIGMAS 
IF(H.LE.20.)GOTO 810 
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INTERPOLATE PR,DR,TR ARRAYS TO GET P,O.T SIGMAS AT HEIGHT H, 
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CALL RTERP(25.,PHI,PR,DR,TR,SPHG, SOHG, STHG) 


GO TO 810 

CONTINUE 
CALL RTERP(H,PHI,PR,DR,TR.SPH, SDH, STH) 
GO TO 820 

LAT-LON INTERPOLATION ON P,D,T SIGMAS AT LOWER HEIGHT 
CALL INTER4( PHi,THET, THA, SP4,504,ST4,PA,DA,TA, 


¢ OPX,.DPY,OTX,OTY) 


LAT-LON INTERPOLATION ON P.D,T SIGMAS AT UPPER HEIGHT 
CALL INTER4{ PHI, THET,IWHB, SP4,SD4,ST4,PB,0B8,.1T8. 


$¢ DPX ,DPY.DTX,.DTY) 
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HEIGHT INTERPOLATION OF SIGMAS 

CALL INTERZ(PA,DA,TA, HA,PB,OB,TB, HB ,SPH,SDH,STH,H) 
IF(PH.LE.0.0.0R.DH.LE.0O.0.0R.TH.LE.O.0)GO TO 825 
IF(H.LE.20.)GOTO 820 


FH = 1. - 0.2*(25. - H) 

SPH = FH*SPHG + (1. - FH)*SPH 
SDH = FH*SDHG + (1. - FH)*SOH 
STH = FHtSTHG + (1. ~ FH)*STH 


HEIGHT DISPLACEMENT BETWEEN PREVIOUS AND CURRENT POSITION 


DZ = Ht - H 

SPHL =SQRT (PLPH*ABS( SPH) ) 

SPHS=SQRT((1.-PLPH)*ABS(SPH) ) 

SDHL =-SQRT(DLPH*ABS( SDH) ) 

SDHS=SORT((1.-DLPH)*ABS(SDH) ) 

STHL=SORT(TLPH*ABS(STH) ) 

STHS=SORT((1.-TLPH)*ABS(STH)) 

SPH = SOQRTCABS(SPH) ) 

SDH = SQRT(ABS(SDH)) 

STH = SORT(ABS(STH)) 

COMPUTES HORIZONTAL DISPLACEMENT DX BETWEEN PREVIOUS AND CURRENT 
POSITION, HORIZONTAL SCALE HL, AND VERTICAL SCALE VL 

COMPUTES PERTURBATION VALUES PRH,ORH,TRH,URH,VRH AND WRH 

CALL PERTRB 

ADDS RANDOM PERTURBATIONS TO PH,DH,TH 

PH = PH*(t. + PRH) 

DHE = DH*(1. + DRH) 

TH = THY( 14. + TRH) 

ADDS RANDOM WINDS TO UH, VH.WH 

UH=UH+URH 

VH=VH+VRH 

WH=WGH+WRH 

SETS PREVIOUS RANDOM PERTURBATION IN P,D,T TO CURRENT 
PERTURBATIONS. FOR NEXT CYCLE 

RP1S= PRHS 

RD1S= DRHS 

RTiS= TRHS 

RP1L=PRHL 

RD 1. =DRHL 

RT 4L=TRHL 

SETS PREVIOUS MAGNITUDES TO CURRENT VALUES, FOR NEXT CYCLE 

SP1S=SPHS 

SD1S= SDHS 

ST1S=STHS 
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SP tL=SPHL 

SO 1L=SDHL 

STAIL=STHL 

SETS PREVIOUS WIND PERTURBATION VALUES TO CURRENT VALUES, 
FOR NEXT CYCLE 

RU1S=URHS 

RV1S=VRHS 

RU 4L=URHL 

RV1L=VRHL 

RW1=WRH 


SETS PREVIOUS WIND PERTURBATION MAGNITUDES TO CURRENT VALUES. 


FOR NEXT CYCLE 
SU 1S=SUHS 
SV1S=SVHS 
SU1L=SUHL 
SV1L=SVHL 
SW1=SWH 
SETS PREVIGUS HEIGHT TO CURRENT HEIGHT. FOR NEXT CYCLE 
Hi = H 
SETS PREVIOUS LATITUDE TO CURRENT LATITUDE, FOR NEXT CYCLE 
PHI 1R=PHIR 
SETS PREVIOUS LONGITUDE TO CURRENT LONGITUDE. FOR NEXT CYCLE 
THETtR=THETR 
SETS NMORE TO COMPUTE MORE DATA ON NEXT CYCLE 
NMORE = 1 
NO MORE DATA IF P, D, OR T LEQ O 
IF (PH*DH*TH.LE.O.) RETURN 
CALL STDATM(H,TS,PS.0S) 
IF ((PS*DS*TS).GT.O.) GO TO 870 
PGHP =O. 
DGHP =O. 
TGHP=0. 
PHP =O. 
DHP=0. 
THP=O. 
GO TO 880 
PGHP= 100. *(PGH-PS)/PS 
DGHP= 100. *(DGH-DS)/DS 
TGHP=100.*(TGH-TS)/TS 
PHP=100.*(PH-PS)/PS 
DHP = 100. *(DH-DS)/DS 
THP=100.*(TH-1S)/TS 
CONVERTS QBO P,D.T TO PERCENT 
PQ=100.*PQ 
DQ=100. *DQ 
TQ=100.*TQ 
CONVERTS RANDOM P,D,T TO PERCENT 
PRH= 100. *PRH 
DRH= 100. *DRH 
TRH= 100. *TRH 
PRHS=100. *PRHS 
DRHS= 100. *DRHS 
TRHS= 100. *TRHS 
PRHL=100. *PRHL 
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579 
980 
581 
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DRHL = 100. *DRHL 
TRHI_.= 100. *TRHL 
SPHS = 100.*SPHS 
SDHS = 100.*SDHS 
STHS = 100.*STHS 
SPHL = 100. *SPHL 
SDOHL = 100. *SOHL 
STHL = 100.*STHL 
CONVERTS WIND SHEAR TO M/S/KM 
DUH = DUH * 1000. 
NVH = DVH * 1000. 
FQA=PQA* 100. 
NQA=DQA* 100. 
[QA=TQAt 100. 
SPH=SPH?* 100. 
SDH=SDH?t 100. 
STH=STH?* 100. 
PSH=PSH* 100. 
DSH=DSH* 100. 
TSH=TSH* 100. 

IF (NPOP .NE.O) GO TO 920 
UPRE=UGH 

VPRE =VGH 

DUPRE =DUH/ 1000. 
DVPRE =DVH/ 1000. 
RETURN 

IF (IOPP.NE.O) 


sf WRITE(IOPP,951)H,PHI.THET ,DGHP, TGH, UGH, VGH, WGH, SOHL ,STHL, 


& SUHL,SVHL 
FORMAT(FS5.1,7F7.2,4F6.2) 


WRITE(6,9C0) H. PHI, THET,PGH,DOGH, TGH, UGH, CHAR(IWSYM), 


1 VGH, PH,DH, TH,UH, CHAR( IWSYM),VH,DUH., 


$¢ DVH,SWH,IET,PGHP ,DGHP , TGHP , WGH, PHP .DHP , THP ,WH,PSH,OSH, TSH, 


$ SPU,SPV,PQ,00,TQ,U0, 


$ VQ,.PQA,DQA,TQA,UA, VA,PRHS,DRHS, TRHS ,URHS, VRHS, SPHS, SOHS. STHS, 
1SUHS .SVHS,PRHL ,ORHL, TRHL,URHL, VRHL ,SPHL ,SOHL,STHL, SUHL.SVHL, 


2PRH, DRH, TRH, URH, VRH, SPH, SDH, STH, SUH, SVH 


FORMAT( 1X, F6.2.2F7.2,2(2E9.3,2F6.0,A1,F5.0),2F5.1,23X,F6.2/1X, 


1 15,14X,2(FB.1,'°%’),FE.1,'°%' ,E10.2, 1X, 
& 2(FB.1,'%').FEG.1,'%', FIO.2, 11%, 
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23F5.1,2F5.0,’ SP’/102X,3F5.1,2F5.0, ° QBO0’/102X,3F5.1,2F5.0,’ MAG’/SCIM 


3 102X,.3F5.1,2F5.0,’ RANS’/102X,3F5.1,2F5.0, ° 


4102X,3F5.1,2F5.0,’ RANL’,/ 
5102X%,3F5.1,2F5.0,’ SIGL’,/ 
6102X%,3F5.1.2F5.0,’ RANT’./ 
7102X,3F5.1,2F5.0,’ SIGT’,/) 
RETURN 

END 


SUBROUTINE SELEC4 
INTEGER IOTEM2 
COMMON/C4/XL(16),YL(€16),NP, DUMMY (2499 ) 


SUBROUTINE TO SELECT POINTS FOR INTERPOLATION 


SIGS’/ 
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COMMON/ IOTEMP/IOTEM1, IOTEM2,DUMMY2(61) 
COMMON /POINT/ IPT(16.5),LL(16),0xXY(16,2) 
COMMON /ORDER/ IPTN( 16,5), IREAD(65,3) 


DIMENSION 10(4),1IL(2),utL(2),LIML(51),L1IMU(51) 


DATA CLIML/15,14,19,12,11,10,9,8,7,6.5,4,9,2,239*1,2,3,4,5,6,7,8.3, 
110,11,12,13, 14, 15/ 
DATA LCIMU/33,34,35,36,37,38,39,40,41,42.43,44,45,46,23*47,46,45, 
144,43,42,41,40,39,38,37.36,35,34,33/ 


DATA PI/3.14159/ 
INITIALIZE 


PI4=PI/4. 
DEGRAD=PI/180. 
DO i 1=1,16 
DO 1 J=1,5 
IPT(I,J)=0 


MAJOR LOOP FOR POINTS 
DO 100 I1=1,NP 


LA=ABS(XL(II))*10.+.5 
LO=YL(IIT)*10.+.5 

IF(LO.LT.O) LO = LO + 3600 

LL( II )=LA* 10000+LO 

IF (XL(IIT).LT.O.) LECII)=-LtLL(TT) 


IF (XL(II)-15.1) 15,30,30 
IF (XL(II)) 50,40,40 


NMC GRID 


IPTC(II,5)=2 
YEL = YL(1I) 
IF(YEL.LT.O.)YEL = YEL + 360. 
EL = (350-YEL)*DEGRAD 
PHI=XL(11)*DEGRAD 
R=31.204359052+(SIN(PI4-PHI/2. )/COS(PI4-PHI/2. )) 
XX=R*COS(EL)+24. 
YY=R*SIN(EL)+26. 
I=XX 
J=YY 
OX=Xx-I 
DY=YY-J 
DxY(1I,1)=Dx 
DOXY(II,2)=DY 
IF (XL(I1I).GT.17.18) GO TO 31 
IF ((JU.LT.1).0R.(JU.GT.51)) GO TO 70 
IF (C(1.LT.LIML(JU)).OR.(1.GT.LIMU(J))) GO TO 70 
I1C(1)=I*1000+J 
IF ((ABS(DX).GT..1).0R.(ABS(DY).GT..1)) GO TO 32 
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34 


35 


38 


40 


45 


50 


IP=1 

GO TO 35 

CONTINUE 

IF (€XLCII).GT.17.18) GO TO 34 


IF ((C1.GT.(LIMU(JU)-1)). AND. ((U.GE.15). AND. (U.LE.37))) 


.OR.(JU.GT.50)) GO TO 70 


IF ((141.GT.LIMU(J+1)).0R.(1.LT.LIML(JU+1))) GO TO 80 
IF ((1.EQ.LIMU(J)).OR.(1I.EQ.LIML(J))) GO TO 80 


IP=4 

I1C(2)=(14+1)* 1000+J 
IC(3)=1* 1000+U+1 
10(4)=(1+1)* 1000+U+1 
CONT INUE 

REWIND IOTEM2 

DO 38 IPG=1,1977 
READ(IOTEM2) TJ 

DO 38 K=1,IP 

IFC IC(K).EQ. IU) IPTC ITI, K)=IPG 
GO TO 100 


EQUATORIAL GRID 


IPT(II,5)=1 
L1=XL(I1) 
L2=YL(II1) 

IF(L2.LT.O)L2 = L2 + 360 
IL(1)=L1/5 

ILC2})=ILC1)+1 
JL(1)=(02/5)+1 
JL(2)=JUL(1)-1 
DO 45 K1=1,2 
DO 45 K2=1,2 

OLON = ABS(YL(II) - 5.*JUL(K2)) 


IF (DLON .GT. 80.) DLON = ABS(360. 


- DLON) 


IF (CABS(XL(II)-5S.*IL(K1)).GI.0.1).0R.(DLON GT.Q.1))GOTO 45 


IF (JL(K2).EQ.72) JL(K2)=0 
IPT(II, 1) =UL(K2)*44+IL(K1 +1 
GO TO 100 

CONTINUE 

IF (JL(1).€Q.72) JL(1)=0 
IPTCIE, 1)=JLC1)*44+ IL (4 +1 
IPT(II,2)=JL(2)*4+I0€1)+1 
IPT(II,3}=UL(1)*4+IL(2)+1 
IPT(I1T, 4) =JL(2)74+IL(2)+1 
GO TQ 100 


SOUTHERN HEMISPHERE 


IPT(I1,5)=3 
Li=XL(II) 
L2=YL( II) 
IF(L2.LT.O)L2 = L2 + 360 
IF (ABS(XL(TI)).LT.85.0) GO TO 51 
IPT(IIE,1)=1 
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IF (ABS(XLCII)4+90.).1.7.0.11) GO TO 
CONT INUE 

IL(1)=(L1/5)-1 

JL(1)=(L2/5)+1 

IL(2)=IL( 1) +1 

JL(2)=UL(1)-1 

DO 52 K1i=1,2 

DO 52 K2=1,2 

DLON = ABS(YL(II) - 5. tJUL(K2)) 


IF(DLON .GT. 180.) DLON = aBs(360. 
IF (CABSEXLCII)-5.*IL(K1)).GT.0.1) 


IF (JL(K2).EQ.72) JUL(K2)=0 
IPT(II, #)=ULEK2)*17-1L(K1)+1 
IF(IL(K1).NE.0}) GO TO 100 
IPTC II, 1) =JL(K2)*41t1 
IPT({II,5)=1 

GO TO 100 

CONT INUE 

IF (JL(1).EQ.72) JL(1)=0 
IF (IPT(II,1).€0.1) GO TO 54 
IPT(II,1)=UL(1)*17-IL01)+1 
IPT({II,2)=JUJL(2)*17-IL( 1)+1 
IF (IL(2)) 55,53,55 
IPT(IIT,3)=Jt (1) *4+1 
IPT(II,4)=JL(2)*4+1 
IPT(II,5)=1133 

GO TO 100 
IPT(II,2)=JUL(1)*17-IL(2)+1 
IPTCII,3)}=JUL(2)*17-3IL02)+1 
IPT(II,5)=333 

GO TO 100 

CONT INUE 
IPT(II,39)=JUL(1)*17-1IL(2)+1 
IPT(II,4)=JUL(2)*17-IL(2)+1 
GO TQ 100 


BORDERLINE POINTS 


CONT INUE 

TWO NMC, TWO EQUATORIAL 
IPT(II,5)=2211 

L=YL(1I1) 
IPT(II.1)=((L/5)+2)*4 
IPT(II,2)=IPT(II,1)-4 

IF (L.GE.355) IPT(II,1)=4 


IF (JU.LT. 14) Jat 

IF (JU.GT.51) J=51 

IF (CI.LT.LIML(J)) T=LIML(J) 
IF (I.GT.LIMU(J)) T=LIMU( Ud) 
1c(1)=I1* 100047 


100 


- OLON) 
.OR.(OLON.GT.O.1))GOTO 52 


IF ((JU.LT.15).OR.(U.GT.37)) GO TO 72 


IC(2)=I1* 1000+J+1 
GO TO 75 
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IF ((JU.NE.1).ANO.(JU.NE.51)) GO TO 74 
IF (1.EQ.LIMU(J)) GO TO 73 
IC(2)=(1+1)*1000+uU 

GO TO 76 

1¢(2)=(1-1)*10004+¥U 

GO TO 76 

IF (I.EQ.LIML(JU)) GO TO 75 
1C(2)=LIMU(J+1)* 1000+UJ+ 1 

GO TO 76 

IC(2)=LIML(JU+1)* 1000+UJ+1 


REWIND IOTEM2 
DO 77 IPG=1,1977 
READ(IOTEM2) IJ 


pO 77 K=1,2 
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IFCIC(K).EQ. IU) IPTC1I.K+2)=1PG 
GO TO 100 


CONT INUE 
THREE NMC, ONE EQUATORIAL 
IPT(II,5)=2212 

Ic(2) = 0 

L=YL(II) 
IPT(T1,2)=((L/5)+1)*4 

IF (L.GE.355) IPT(II,2)=4 
IF (1.£€Q.LIML(J)) GO TO 84 
IF (J.GT.37) GO TO 82 
IC(1)=I*1000+d 

IC(3)=I1* 1000+J+1 
1C(4)=(14+1)}* 1000+UJ+1 

GO TO B88 

IC€4)=(1+1)* 1000+¥U 
IC(3)=1* 10004J 
IC(4)=1*1000+U+1 

GO TO 88 

IF (J.GT.37) GO TO 86 
Ic(t1)=(1-1)*1000+uU+1 
IC(3)=I*1000+UJ+1 
1C(4)=1+*1000+J 

GO TO 88 

IC( t)=(1+1)*1000+U+1 
IC(3)=(1+1t)*1000+uU 
IC(4)=1* 1000+J 


REWIND IOTEM2 

DO 839 IPG=1, 1977 

READCIOTEM2) IJ 

DO 89 K=1,4 

IF(IC(K).EQ.0) GO TO 89 
IF(IC(K).€Q.IJU) IPT(II,K)=IPG 
CONTINUE 


CONT INUE 
DO 150 I=1,16 
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... 1TOPR=414 RANDOM OUTPUT 


... IOPQ=f QBO OUTPUT 
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DO 150 J=1,5 


150 IPTN(I,JU)=IPT(1I,J) 


CALL SORT4(NP) 

RETURN 

END 

SUBROUTINE SETUP 
CHARACTER*2 IC 
COMMON /COTRAN/NDATA(19),MI,IH,IEX,1x( 10) 

COMMON /COTRAN1/IC 

COMMON /IPRTP/IPRT,NLIMIT 
DIMENSION IP(5),ID(5),I1T(5), 1DAY(12),xXWR(29) 

COMMON /IOTEMP/IOTEM1, IGTEM2,1US,0D,XMUD,PHI1,PHI, 
NSAME ,RPIL,RDOIL,RTIL,SPIL, SOIL, STIL, RUIL,RVIL,SUIL,SVIL, 
$¢ MN, IOD, IYR, H1, PHIIR, THETAIR,DUMS(21),RP1S,RD1S 
1,RT1S,RUIS,RVIS,SP1S,5015,5T1S,5U1S,5V1S,UDS1,VDS1, 
2UDL1,VDL1,U0S2,VDS2,UDL2,VDL2 

COMMON /PDTCOM/1U4,MONTH, IOPR, 

PSP(15,19,18),DSP(15,19,18),TSP(15,19,18),USP(15.19,18), 
VSP(15,19,18), 
PG(21,19),.0G(21,19),TG(21,19),UG(21, 19), 
PAQ(17,5),0AQ(17,5),TAQ(17,5),UAQ( 17,5), VA0(C 17,5), 
PDQ@(17.5),000(17,5),T0Q(17,5),UD0(17,5),VvD0( 17,5), 
— PR(29,19),0R(29,19),TR( 29, 19),UR(29,19),VR(29, 19), 
* PQ,DO,TQ,UQ, VQO,PQA,DQA, 

. TQA,UA,VA, IOPQ,PLP(25,10),DLP(25, 10), TLP(25, 10) 
1,ULP(25,10),VLP(25, 10) ,UOL(25,10),VDL(25,10) ,UDS(25, 10) 
2,V0S5(25, 10) 

COMMON /CHIC/OUM( 18), I1WSYM,USH, VSH,DUSH,DVSH 

COMMON /RAND/IXX,IY,12 

COMMON /VERT/RW1,SW1,WRH, SWH,WR(29) 

DIMENSION IDUM(9) 

DATA IDAY/0,31,59,90, 120, 151, 181,212,243,273,304 ,334/ 
DATA XWR/.87,.15,.14,.08,.04, .04, .05, .07,.11,.14,.22, 
&® .34,.46,.73,1.30,2.11,3.40,6.30,9.58,10.55,11.51,11.32,11.13, 
&® 10.94,10.75,10.85,11.76,11.86,12.06/ 
po 141 = 1,29 

WR(I) = XWR(T) 

XMUD = O. 

IF (MN.GT.12) GO TO 2 

IDA = IDAY(MN) + IOD 


DD = IDA . 

IF (MOD(IYR,4).EQ.0.AND.MN.GT.2) IDA = IDA + 1 

XMJD = 2439856. + 365. * (IYR - 68.) + IDA + INT(CIYR - 65.) 
$ / 4.) 


...S5SECOND DATA LINE READS, FREE FIELD, THE FOLLOWING DATA: 


IUS = UNIT NUMBER FOR SCIDAT DATA 

IU4 = UNIT FOR 4-D INPUT P,D,T O-25KM DATA 

IOPR = RANDOM OUTPUT OPTION 
IOPR=2 NO RANDOM OUTPUT 
IOPQ = OBO OUTPUT OPTION 
IOPQ=2 NO QBO OUTPUT 
NR1 = STARTING RANDOM NUMBER 


... TOTEM1=UNIT FOR 4-D P, D, T DATA (SCRATCH FILE, DOES NOT NEED TO 


BE ASSIGNED) 
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sists. IOTEM2=UNIT FOR NMC GRID POINTS (SCRATCH FILE, DOES NOT NEED TO SETU 
BE ASSIGNED) SETU 

IPRT = O FOR 4-D DIAGNOSTIC OUTPUT, = 1 FOR NONE SETU 

NLIMIT = MAXIMUM NUMBER OF TEST VIOLATIONS ALLOWED IN 4D SETU 

DATA (12 LE NLIMIT LE 25) SETU 

2 READ(5,*) IUS,1U4,IOPR, IOPQ,NR1,IOTEM1, IOTEM2,IPRT,NLIMIT SETU 
IFCNLIMIT.LT.12) NLIMIT=12 SETU 
IFCNLIMIT .GT.25) NLIMIT=25 SETU 
WRITE(6,9000) IUS,IU4,IOPR, IOPQ,NR1,IOTEM1, IOTEM2 SETU 

¢ ,XMJD SETU 

IF (IOPR.LT.1.OR.IOPR.GT.2) GO TO 666 SETU 

IF (IOPQ.LT.1.0R.IOPQ.GT.2) GO TO 666 SETU 
MONTH=MN SETU 
RPSCALE = 1.0 SETU 

IF (IOPR.EQ.2) GO TO 7 SETU 
IF(NRt . LE. O .OR. NR1 .GE. 30000)STOP ‘FIRST RANDOM NUMBER OUT OFSETU 

¢ RANGE ’ SETU 

IXX = NR1 SETU 

IY = 172 * MOO(IXX, 176) - 35 * €IXxX / 176) SETU 

IZ = 170 * MOD(IXX, 178) - 63 * (IXX / 178) SETU 

IF (IY .LT. O) IY = IY + 30307 SETU 

IF (IZ .LT. O) 12 = 12 + 303923 SETU 

R = RANDOM(L) SETU 

IF(L .£Q. 1)STOP ‘BAD CALL TO RANDOM ON FIRST TRY’ SETU 

Sag Te as THIRD DATA LINE READS FREE FIELD, THE FOLLOWING DATA SETU 
FACH FOR LARGE SCALE (L)} AND SMALL SCALE (5S) COMPONENTS SETU 
RPiL,RP1S = INITIAL RANDOM PRESSURE PERTURBATION, PERCENT SETU 
RDiL,RD1S = INITIAL RANDOM DENSITY PERTURBATION, PERCENT SETU 
RTIL,RT1S = INITIAL RANDOM TEMPERATURE PERTURBATION, PERCENT SETU 
RU1IL,RUIS = INITIAL EASTWARD WIND PERTURBATION, W/S SETU 
RVIL,RV1S = INITIAL NORTHWARD WIND PERTURBATION, W/S SETU 
RPSCALE = RANDOM PERTURBATION SCALE, NOMINAL=1.0 SETU 
MAXIMUM=2.0 , MINIMUM=-0.0 SETU 

READ(5,*)RPIL,RP1S ,ROIL,RDIS,RTIL,RT1IS,RUIL,RUIS, SETU 

& RV1L,RV1S,RPSCALE SETU 
IF((RPSCALE .LT.0.0).0R.(RPSCALE .GT.2.0)) RPSCALE=1.0 SETU 
RW1=0. SETU 
RP1=RPIL+RP1S SETU 
RD1=RD1S+RDIL SETU 
RTI=RTIStRTIL SETU 
RUT=RUIL+RUIS SETU 
RVI=ERVIL+RVIS SETL 
AVOIDS FILE SEARCH IF CURRENT MONTH IS SAME AS PREVIOUS MONTH SETU 
IF(NSAME.GT.O) GO TO 621 SETU 

7 IF (NSAME.EQ.1+) GO TO 621 SETU 
CALL GETNMC SETU 

sei ete LOADS NMC GRIO DATA FROM INPUT UNIT TO SCRATCHFILE UNIT IOTEM2 SETU 
IF (MONTH .LT. 1 .OR. MONTH .GT. 12)GOTO 666 SETU 

13 DO 100 1=1,252 SETU 
CALL RTRANt SETU 
READS ZONAL PRESSURE DATA SETU 
IF(IC .NE. ‘2P’ }GOTO 666 SETU 
IF(MI .NE. MONTH)GO TO 100 SETU 

50 IH=(IH-15)/5 SETU 


*_ es © 8 # 


304 
305 
308 


320 


350 


360 


TENX=10.**1EX 

DO 60 J=1,19 

PG(IH, J) =NDATA(J) *TENX 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONT INUE 

DO 200 I=1,252 

CALL RTRAN1 

READS ZONAL DENSITY DATA 

IF(IC .NE. ’ZD’)GO TO 666 

IF(MI .NE. MONTH)GO TO 200 

IH=(1IH-15)/5 

TENX=10.**IEX 

D0 160 J=1,19 

OG(1IH,J)=NDATA(J)*TENX 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

DO 300 I=1,252 

CALL RTRAN14 

READS ZONAL MEAN TEMPERATURE DATA 

IF(IC .NE. ‘2T’)GO TO 666 

IF(MI .NE. MONTH)GOTO 300 

IH=(IH-15)/5 

TENX=10.**IEX 

DO 260 J=t,19 

TGC IH, J) = NDATA(JU)*TENX 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONT INUE 

pO 305 I=1,252 

CALL RTRANI 

READS ZONAL AVG. ZONAL WIND DATA 

IF(IC .NE. ’ZU’)GO TO 666 

IF(MI .NE. MONTH)GOTO 305 

IH=(IH-15)/5 | 

TENX=10. ** TEX 

DO 304 J=1,19 

UG( IH, J)=NDATA(UJU)*TENX 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

DO 360 I=1,3060 

CALL RTRAN 

READS STATIONARY PERTURBATION DATA FOR PRESSURE (TO BE STORED IN 
PSP ARRAY) 

IF(IC .NE. ‘SP’)GO TO 666 

IF(MI .—£Q. MONTH) GO TO 320 

GO TO 360 

ISH=(IH-15)/5 

K=(NDATA(t)+100)/10 

DO 350 L=1,18 

PSPCISH,.K,L)=NDATA(L+1)/1000. 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONT INUE 

DO 365 I=1,3060 

CALL RTRAN 

READS STATIONARY PERTURBATION DATA FOR DENSITY (TO BE STORED IN 
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SETU 
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65 


66 


362 


364 


365 


367 


369 


370 


372 


374 


375 


377 


379 


380 


DSP ARRAY) 

IF (IC .NE. ‘SD’)GO TO 666 

IF (MI .EQ. MONTH)GOTO 362 
GO TO 365 

ISH = (IH-15)/5 

K = (NDATA(1)+100)/10 

DO 364 L=1,18 
DOSP{ISH,K,L)=NDATA(L+1)/1000. 
CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

00 370 I=1,3060 

CALL RTRAN 


SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 


READS STATIONARY PERTURBATION DATA FOR TEMPERATURE(TO BE STORED INSETU 


TSP ARRAY ) 

sF (IC .NE. ‘ST’)GO TO 666 

IF (MI .EQ. MONTH)GOTO 367 
GO TO 370 

ISH = (IH-15)/5 

K = (NDATA(1)4+100)/10 

DO 369 L=1,18 
TSP(CISH,K,L)=NDATA(L+1)/1000. 
CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

DO 375 1r=1,3060 

CALL RTRAN 


READS STATIONARY PERTURBATION DATA FOR ZONAL WIND(TO BE STORED IN 


USP ARRAY ) 

IF (IC .NE. ‘SU’)GO TO 666 
IF (MI .EQ. MONTH)GOTO 372 
GO TO 375 

ISH = (IH-15)/5 

K = (NDATA(1)+100)/10 

DO 374 L=1, 18 
USP(ISH,K,L)=NDATA(L+1)/10. 
CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

DO 380 I1=1,3060 

CALL RTRAN 


READS STATIONARY PERTURBATION DATA FOR MERIDIONAL WINO(TO BE 


STORED IN VSP ARRAY) 

IF (IC .NE. ‘SV’)GO TO 666 

IF (M1 .EQ. MONTH)GOTO 377 

GO TO 380 

ISH = (IH-15)/5 

K = (NDATA(1)+100)/10 

DO 379 L=1,18 

VSPCISH,K,L)=NDATA(L+1)/10. 

CONVERSION TO REAL AND STORAGE IN ARRAY COMPLETE 
CONTINUE 

ZERO STATIONARY PERTURBATIONS AT LATITUDE +90 AND -90 
DO 382 1=1,15 

DO 382 L=1.18 

DO 382 K=1,19,18 

PSP(I,K,L)=0. 


SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 


158 
159 
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382 


410 


420 
430 


435 


A440 
445 


450 


455 
460 


465 


470 
475 


DSP(1,K,L) 
TSP(I,K.L) 
USP(I,K,L) 
VSP(I,K,L) 
CONT INUE 
DO 430 I=1,300 

READS RANDOM PERTURBATIONS IN PRESSURE 

CALL RTRAN1! 

IF(IC .NE. ‘RP‘)GO TO 666 

IF(MI .EQ. MONTH)GO TO 410 

GO TO 430 

IHR = (IH+5)/5 

IF(CIH .GT. 120) IHR=194+IH/20 

DO 420 K=1,19 
PR(IHR,K)=(NDATA(K)*RPSCALE/ 1000 .0)**2 
CONTINUE 

DO 445 I=1,300 

READS RANDOM PERTURBATIONS IN DENSITY 

CALL RTRANT 

IF(IC .NE. ’RD’)GO TO 666 

IF(MI .EQ. MONTH)GO TO 435 

GO TO 445 

IHR = (IH+5)/5 

IF(IH .GT. 120)IHR=194+IH/20 

DO 440 K=1,19 

DRC IHR. K)=(NDATA(K)*RPSCALE/ 1000 .0)**2 

CONT INUE 

DO 460 I=1,300 

READS RANDOM PERTURBATIONS IN TEMPERATURE 
CALL RTRAN1 

IF(IC .NE. ‘RT’)GO TO 666 

IF(MI .EQ. MONTH)GO TO 450 

GO TO 460 

IHR = (IH+5)/5 

IF(IH .GT. 120)IHR=19+IH/20 

DO 455 K=1,19 

TRC IHR, K)=(NDATA(K)*RPSCALE/1000.0)**2 

CONT INUE 

DO 475 I=1,348 

READS RANDOM PERTURBATIONS IN ZONAL WIND 
CALL RTRAN{ 

IF(IC .NE. ‘RU’)GO FO 666 

IF(MI .—EQ. MONTH)GO TO 465 

GO TO 475 

IHR = (IH+5)/5 .- 

IF(IH .GT. 120) IHR=19+IH/20 

DO 470 K=1,19 
URCIHR.K)=(NDATA(K)*RPSCALE/10.0)**2 

CONT INUE 

DO 490 1I=1,348 

READS RANDOM PERTURBATIONS IN MERIDIONAL WIND 
CALL RTRAN1 

IF(IC .NE. ‘’RV’)GO TO 666 

IF(MI .EQ. MONTH)GO TO 480 
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490 
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810 


820 


830 


840 


850 
855 


860 
B65 


875 
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GO TO 490 

IHR = (IH+5)/5 

IF(IH .GT. 120) ITHR=19+IH/20 

DO 485 K=1,19 

VRC IHR, K)=(NDATA(K)*RPSCALE/10.0)**2 
CONTINUE 

RANDOM SIGMAS ARE ZEROED IF IOQPR=2 
IF(TOPR .EQ. 1)GO TO 500 

DO 495 [=1,29 

DO 495 J=1,19 


DO 840 I=1,25 
CALL RTRAN2 
DO 810 K=1,5 
IP(K)=NDATA(2+K) 
ID(K)=NDATA(7+K) 
IT(K)=NDATA( t2+K) 
IF(IH.GT.90) IH=70+(1H/4) 
IH=t+(IH/5) 
IF(IC .NE. ‘P ‘ .OR. IH .NE. I)GOTO 666 
DO 830 J=1,5 
PLP(I,J+5)=IP(JU)/1000. 
PLP(I,6-J)=IP(JU)/1000. 
DLP(1I,J+5)=1D(JU)/ 1000. 
DLP(I,6-J)=ID(JU)/ 1000. 
TLP( IT, J+5)=I1T(JU)/1000. 
TLP(I,6-JU)-=I1T(J)/1000. 
CONTINUE 
DO 865 1=1,25 
CALL RTRAN2 
DO 850 K=1,5 
IP(K)=NDATA(2+K) 
IO(K)=NDATA(7+K) 
IF(IH.GT.90) IH=70+(IH/4) 
ITH=1+(1H/5) 
IFCI.NE.IH.OR.IC.NE.’PW') GO TO G66 
DO 860 J=1,5 
ULP(I,J+5)=IP(U)/ 1000. 
ULP(1,6-JU)=IP(U)/1000. 
VLP(I,J+5)=ID(JU)/1000. 
VLP(1,6-JU)=I10(JU)/1000. 
CONT INUE 
DO 888 171,25 
CALL RTRAN2 
DO 875 K=1,5 
IP(K)=NDATA(2+K) 
ID(K )=NDATA(7+K) 
IF(IH.GT.90) IH=70+(IH/4) 
IH=1+(IH/5) 
IF(IH.NE.1I.0R.IC.NE.’CS’) GO TO 666 
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527 
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DO 885 J=1,5 
UDS(1I,J+5)=(IP(JU)/ 1000. 
UDS(1,6-J)=(IP(JU)/ 1000. 
VOS(1I,u+5)=(1D0(U)/ 1000. 
VOS(1,6-J)=(ID(JU)/ 1000. 
CONTINUE 
DO 898 1=1,25 
CALL RTRAN2 
DO 892 K=1,5 
IP(K )=NDATA(2+K) 
ID(K)=NDATA(7+K) 
IFCIH.GT.90) IH= 70+(IH/4) 
IH=1+(1IH/5) 
IF(IH.NE.I.OR.IC.NE.’CL’) GO TO 666 
DO 896 J=1,5 
UDL(I,J+5)=(IP(JU)/1000. ) 
UDL(I,6-U)=(IP(U)/ 1000. ) 

«) 

.) 


i 


VOL(1I,J+5)=(10(U)/ 1000 
VDL(1I,6-JU)=(1D0(U)/ 1000 
CONT INUE 

IF(IOPR .EQ. 1)GO TO 910 
DO 905 1=1,25 

DO 905 J=1, tO 
PLP(I,JU)=0O. 
DLP(I,J)=0. 
TLP(I.JU)=O. 
ULP(I,JU)=0O. 
VLP(1,JU)=0. 
UDS(I,uU)=0. 
UDL(I,JU)=0. 
VDS(I,JU)=0. 
VOL(1I,J)=0. 

CONTINUE 

DO 530 1=1,16 

CALL RTRAN2 

IH=MI 

IF (IC.NE.’QP’) GO TO 666 
IH = (IH-5)/5 

DO 530 J=1,5 


CONVERT FROM INTEGER PER MIL - QBO PRESSURE AMPLITUDE 


PAQ(IH,J) = IX(2*U-1)/1000. 

QBO PRESSURE PHASE (DAYS PAST JAN O, 1966) 
POQ(IH,J) = IX(2*JU)*1. 

DO 531 I = 1,5 

PAQ(1,I) = O. 

CALL PHASE(PDQ(2.1),15..PD0(3,1),20.,PDQ(1,1),10.) 
DO 540 I=1,16 

CALL RTRAN2 

IH=MI 

IF (IC.NE.’QD’) GO TO 666 

IH=(IH-5)/5 

DO 540 J=1,5 


C...CONVERT FROM INTEGER PER MIL - QBO DENSITY AMPLITUDE 


DAQ(IH,J) = IX(€2*JU-1)/ 1000. 
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69 


70 


541 


947 


567 


QBO DENSITY PHASE (DAYS PAST JAN 0, 1966) 
DDQ(C IH, J) =1IX(2*dU) 4. 

DO 541 I = 1,5 

DAQ( 1,1) = O. 

CALL PHASE(DDQ(2,1),15.,D0Q0€3,1),20. .D0Q(1.1), 10.) 
DO 550 1=1,16 

CALL RTRAN2 

IH=MI 

IF (IC.NE.‘QT’) GO TO 666 

IH = (IH- 5)/5 

DO 550 J=1.5 


.CONVERTS FROM INTEGER PER MIL - QBO TEMPERATURE AMPLITUDE 


TAQ(IH.JU) = IX(€2*Jd-1)/1000. 

Q@BO TEMPERATURE PHASE 

TOQ(IJH,JU) = IX(2*JU)*1. 

DO 551 I = 1,5 

TAQ(1,1) = O. 

CALL PHASE(TDQ(2.1).15.,.,1T0Q(3,1),20. ,.TOQ(1,1),10. ) 
DO 560 I=1,16 

CALL RTRAN2 

IH=MI 

IF (IC.NE.’QU’) GO TO 666 

IH=(IH- 5)/5 

DO 560 J=1,5 

EASTWARD WIND QBO AMPLITUDE - CONVERTED TO M/S 
UAO(IH,J) = IX(2 * J - 1) / 10. 

EASTWARD WIND Q8O PHASE (DAYS PAST JAN O, 1966) 
UDQ(C IH, J) =IX(2*U)*1. 

DO 564 I = 1,5 

UAQ(1,I) = O. 

CALL PHASE(UDQ(2,1),15.,U0Q(3,1)},.20. ,.uUDQ( 1.1), 10. ) 
DO 570 I=1,16 

CALL RTRAN2 

IH=MI 

IF (IC.NE.‘QV’) GO TO 666 

IH=(IH- 5)/5 

DO 570 J=1,5 

NORTHWARD WIND QBO AMPLITUDE - CONVERTED TO M/S 
VAQ(IH,JUJ) = IX(2 * J - 1) / 10. 

NORTHWARD WIND OBO PHASE (OAYS PAST JAN O, 1966) 
VDO(IH, J)=IX(2*u)*1. 

DO 571 I = 1,5 

VAQ(1,1) = O. 

CALL PHASE(VDQ(2,1),15.,V00(3,1),.20. ,.v0Q(1,1),10. ) 
IF({IOPQ .EQ. 1)GOTO 611 

ZEROS QBO PARAMETERS IF IOPQ=2 

00 610 I=1,17 

DO 610 J=i1, 
PAQ(I,J) 
DAQ(I,uU) 
TAQ( I,J) 
PDO(I,uU) 
DDO(1I,J) 
TOOCI.J) 


Q000008 


fn | |) | 
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407 
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Cc 


UAQ( I,J) =O. 
UDQ(I,J)=0. 
VAQ(I,JU)=0. 
VOQ(I,JU)=0. 


610 CONTINUE 
Cas 


611 


62 1 


REWINDS FILE UNIT IUS 
REWIND IUS 


R=H 1 

IF(H1.LT.25.) R=25. 

CALL RTERP(R ,PHI!,PR,DR,TR.SP1,.SD1,.ST1) 
CALL INTR25(PLP,DLP.H1,PHI1,PLP1,DLP1) 
CALL INTR25(TLP,DLP,H1,PHI1,TLP1.R) 
CALL INTRW(WR,H1,SW1) 

SP 1L=SQRT(PLP1*ABS(SP1))* 100. 

SP 1S=SORT((1.-PLP1)*ABS(SP1))*100. 
SD1L=SQRT(OLP1*ABS(SD1))* 100. 

$0 1S=SQRT((1.-DLP1)*ABS(SD1))* 100. 
ST1IL=SQRT(TLP1*ABS(ST1))* 100. 
ST1S=SORT((1.-TLP1)*ABS(ST1)})*100. 
CALL INTRUV(UR,VR,H1,PHI1,SU1,SV1) 
CALL INTR25(ULP,VLP,H1,PHI1,ULP1,VLP1) 
SUIL=SORT(ULP1*ABS(SU1) ) 
SU1S=SQRT((1.-ULP1)*ABS(SU1)) 
SViIL=SQRT(VLP1*ABS(SV1)) 
SV1S=SQRT((1.-VLP1)*ABS(SV1)) 

CALL INTR25(UDL,.VDL,H1,PHI1,U0L1,VDL1) 
CALL INTR25(UDS,VDS,H1,PHI1,UDS1,V0S1) 
UDL 1=UDL1* 100. 

VDL 1=VDL1* 100. 

UDS1=UDS1t* 100. 

VDS 1=V0S1*100. 


WRITE(6,9001) RPIL,RDIL,RT1IL,SPIL,SDIL.STIL,RUIL .RVIL,SUIL,SVIL, 


1‘ LARGE’ 
WRITE(6,9001)RP1S,RDIS,RT1S,SP1S,S5015.S5T1S,RUIS,RV1S 
1,5U1S,S5V1S, “SMALL’ 
WRITE(6,9002)U0L1,VDL1,UDS1,.VDS1 
WRITE(6,9870)RPSCALE ,IPRT,NLIMIT 
WRITE(6,9003) 

RP1L=RP1L/100. 

RDIL=RDIL/ 100. 

RTIL=RT1L/100. 

SP1L=(SP1L/100. ) 
SO1L=(SD1L/100. ) 
STiIL=(ST1L/100..) 

RP1S5=RP1S/100. 

RD1S=RD1S/100. 

RT1S=RT1S/100. 

SP 1S$=SP1S/100. 

$D1S=SD1$/100. 

S$T1S=S$T1$/100. 

UDL1=UDL 1/100. 

VDL1=VDL1/100. 

UO0S$1=UDS1/100. 
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SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 
SETU 


428 
429 
430 
43 
432 
433 
A34 
435 
436 
A437 
438 
439 
AAO 
aAal 
442 
443 
444 
445 
A4G 
A447 
448 
449 
450 
AS 1 
452 
453 
454 
A55 
A456 
AST 
A458 
459 
A460 
461 
462 
463 
AG4 
465 
466 
467 
468 
469 
470 
474 
472 
473 
ATA 
475 
A476 
A77 
478 
479 
480 
4B 1 


71 


72 


666 


700 


630 FORMAT(27X,’UNPERTURBED (MONTHLY MEAN)’,11X, “MEAN PLUS PERTURBATIOSETU 


INS’ ,9X,’ THERMAL’ ,/,23X,2(34(’-'),2X),93X, ‘WIND’ ,6X, ‘PERTURBATION VASETU 
2LUES’,/.‘ HEIGHT LAT WEST PRES. DENS. TEMP MEAN SETU 

3 PRES. DENS. TEMP TOTAL SHEAR’ ./,2X, ‘(KM)’, 11X, ’LOSETU 

4N’ ,4x%,°(NT/ (KG/ (DEG WIND (M/S) (NT/ {KG/ (DEG SETU 
5WIND (M/S) (M/S/KM) %,28(‘-’),/,.° TIME (DEG) (DEG)’,2¢’ M*t*SETU 

62) M*¥*3) KEL= “.100*=")),2x,6C*=")."7 . P D T U V SETU 

7 S$1GW’/’ (SEC)’,35X,’VIN) E~-W N-S’,20X,’VIN) E-W N-S  E-W SETU 
BN-S (%) (%) (%) M/S M/S M/S'/45X,2(05X, ‘VERT.’ ,26X)) SETU 
9000 FORMAT(’ SCIDAT INPUT UNIT = ‘,12,743,’4-D0 INPUT UNIT = ’,12,/ SETU 
& |’ RANDOM OPTION = ’,12,743,’QB0 OPTION = ’,12,/, SETU 

& * FIRST RANDOM NUMBER = ’,15,/, SETU 

& ’ 4-D P,D,T DATA SCRATCH UNIT = ’,12,743, SETU 

& 'NMC GRID POINTS SCRATCH UNIT = ‘’,12,1783,’JULIAN DATE = ‘’, SETU 
SF9.1,/) SETU 
9001 FORMAT(’ INITIAL P,D,T = ’',3(F6.2,’ % ‘),T6O, ‘SIGMA P,D,T = ’, SETU 
19(FG.2,° % ‘Y,/,° INITIAL ULV = ’,2(F7.2,' M/S ‘),T60, ‘SIGMA SETU 
2U,V = °,2(F7.2,'° M/S ‘yy, 7X,A5,1X,‘SCALE’/) SETU 
9003 FORMAT(//’ *+ PERCENT DEVIATIONS FROM 1976 US STANDARD ' SETU 
1 ,’ATMOSPHERE APPEAR BELOW PRESSURE, DENSITY AND TEMPERATURE ’, SETU 

2 ‘VALUES **/‘//) SETU 
9002 FORMAT(’ INITIAL UDL,.VDL = ‘,2(F6.2,' % ‘), SETU 
iTGO, ‘INITIAL UDS,VDS = ‘,2(FG6.2,’ % ’)) SETU 
9870 FORMAT(/’ RANDOM PERTURBATION SCALING FACTOR = ’,F7.3,T60, SETU 
&® ‘4-0 DIAGNOSTICS PRINT OPTION = ‘,12,/, SETU 

& ° MAX. NUMBER OF 4-D TEST VIOLATIONS = ‘,14) SETU 

END SETU 

SUBROUTINE SLV ( DEN , ALT , XLAT , DAY ) SLY 
SLV 


Cr+ 
C*+ 
C* * 
Cx 
C*r 
Cr 
C++ 
ctr 
C** 
cr 
cC+r 
C+ 
c*, 
cyr 
C+. 
Cyr * 
Cre 


VDS1=VDS1/100. SETU 
WRITE(6,630) SETU 
RETURN SETU 
WRITE(6, 700) IUS,IOPR, IOPQ,NR1,IOTEM1, IOTEM2, SETU 
$¢MONTH,IC,MI,IH, IX, TEX, IP,ID.IT,SO1,NDATA SETU 
FORMAT(’ ERROR IN SETUP INPUT’ ,/,315,110,319,1X,A2,13,14./,4114, SETU 
$/,1514,/,F10.1,/,1916) SETU 
STOP SETU 


i i Na A a a aE a a a i a a a En oe in i ne OA 


‘SLV’ COMPUTES THE SEASONAL-LATITUDINAL VARIATION OF DENSITY IN THE*SLV 
LOWER THERMOSPHERE IN ACCORDANCE WITH L. JACCHIA, SAO 332, 1971. **SLV 
THIS AFFECTS THE DENSITIES BETWEEN 90 AND 170 KM. ‘SLV’ NEED NOT **SLV 
CALLED FOR DENSITIES ABOVE 170 KM, BECAUSE NO EFFECT IS OBSERVED. **SLV 

TSUN 
THE VARIATION SHOULD BE COMPUTED AFTER THE CALCULATION OF DENSITY **SLV 
DUE TO TEMPERATURE VARIATIONS AND THE DENSITY(DEN) MUST BE IN THE **SLV 
FORM OF A BASE 10 LOG. NO ADJUSTMENTS ARE MADE TO THE TEMPERATURE**SLYV 


OR CONSTITUENT NUMBER DENSITIES IN THE REGION AFFECTED BY THIS 7 ESL 
VARIATION. *ASLYV 
Te SkEyY 

DEN = DENSITY (LOGI1O) Re SY 

ALT = ALTITUDE (KM) **SLV 

XLAT = LATITUDE (RAD) **SLV 

DAY = DAY NUMBER TT SLV 

saibin f 


482 
483 
484 
485 
486 
487 
488 
489 
4390 
49} 
492 
493 
494 
495 
496 
497 
498 
499 
500 
501 
502 
503 
504 
905 
506 
507 
508 
309 
510 
511 
So t2 
513 
514 
5915 
516 


OMA DUA WAN a 


10 
11 
12 
Ls 
14 
aes 


17 
18 
19 


CR RT RE EEE ETE ET OEE OE EERIE ESL Y 


C** INITIALIZE DENSTTY (DEN) = 0.0 
c 
DEN = 0.0 
C 
C** CHECK IF ALTITUDE EXCEEDS 170 KM 
C 
IF ( ALT. GT. 170. ) RETURN 
C 
C** COMPUTE DENSITY CHANGE IN LOWER THERMOSPHERE 
C 
Z = ALT - 90. 
X = -0.00139 * Z * Z 
Y = 0.0172 * DAY + 1.72 
P = SIN (Y) 
SP = ( SIN (XLAT) ) **2 
S = 0.014 * Z * EXP (X) 
D S * P * SP 
Cc 
C** CHECK TO COMPUTE ABSOLUTE VALUE OF ‘XLAT’ 
C 
IF ( XLAT. LT. O. ) D = -D 
DEN = D 
RETURN 
END 
SUBROUTINE SLVH ( DEN , DENHE XLAT SDA ) 


SLV 
SLV 
SLV 
SLV 
SLY 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLYV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLV 
SLVH 
SLVH 


CII IRIN RRR ROR ROR IOI ROHR IO III IO RT SEE EEE RHEE HTT ESSE EF Y SL VE 


C* * 
C** NUMBER DENSITY ACCORDING TO 
C** CORRECTION IS NOT IMPORTANT 
C** 


C** DEN = 
c** DENHE = 
c** XLAT = 
c* SDA = 


‘SLVH’ COMPUTES THE SEASONAL-LATITUDINAL VARAITION OF THE HELIUM 


L. JACCHIA, SAO 332, 1971. THIS 


BELOW ABOUT SOO KM. 


DENSITY (LOG10O) 
HELIUM NUMBER DENSITY (LOG10) 
LATITUDE (RAD) 
SOLAR DECLINATION ANGLE (RAD) 


**SCVi 
**SLVH 
**SLVH 
**#SLVH 
**SLVH 
~* SUV 
**SLVH 
**SLVH 


CERES IORI III RII OORT RE SHEE ERE TEEN ESE RSH ESL VH 


©.40909079 ) ) 


B 


DO = 10. ** DENHE 

A = ABS ( 0.65 * ( SDA / 

B = 0.5 * XLAT 
C 
C** CHECK TO COMPUTE ABSOLUTE VALUE OF ‘B’ 
Cc 

IF ( SDA. LT. O. ) Be - 
C 
C** COMPUTE X, Y, DHE AND DENHE 
C 

X = 0.7854 - B 

Y = ( SIN (X) ) ** 3 

DHE= A * ( Y - 0.353956 ) 


SLVH 
SLVH 
SLVH 
SULVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 


73 


74 


C 


DENHE = DENHE + DOHE 


C** COMPUTE HELIUM NUMBER DENSITY CHANGE 


C 


oi ee & ee 


on) 


6 oh i Oe 
DEL= Dt - 
RHO= 10. 


** DENHE 
DO 


++ 


DEN 


DRHO = ( 6.646E-24 ) * DEL 
RHO = RHO + DRHO 


DEN 


RETURN 
END 


ALOG1O (RHO) 


SUBROUTINE SORT4(NP) 


SORTS POINTS FOR SEQUENTIAL FILE READING 


ASSIGNS POINT NUMBERS BY 


COMMON /ORDER/ IPT (16,5), IREAD(65,3) 


DO 1 1=1,65 
DO 1 J=1,3 


IREAD(I,JU)=0 


DO 9 I=1t,NP 


IFC IPT(CI,S). 
IFCIPT(I,5). 
IF(IPT(1,5). 
IF(IPT(1.5). 
IF(IPT(I,5). 
IFCIPT(I,5S). 
.EQ. 


IFC IPT(I,5S) 


.1) GO TO 10 
.1) GO TO 9 
.2) GO TO 2 
.3) GO To 4 
.1133)GO TO 6 
.2211) GO TO 7 


2212)GO TO 8 


IF (IPT(I,5).£0.333) GO TO 4 


GO TO 10 
DO 3 J=1,4 


IFCIPT(I,J).LT.1) GO TO 3 


IPT(I,J)=IPTCI, 


CONTINUE 
GO TO 9 
DO 5 J=1,4 


IFC IPT(I,dJ). 


LT. 
,J)+2265 


IPT(I,J)=IPT(I 


CONTINUE 
GO TO 9 
IFC IPT(I, 1) 
IFC IPT(I,2) 
GO TO 9 
IFC IPT(1,3) 
IFC IPT(1I,4) 
GO TO 9 
IFCIPT(I,1) 


IF(IPTCI,3). 
IFC IPT(I,4). 


CONTINUE 


GT. 
.GT 


.GT 
.GT. 


.GT 


GT 
GT 


J)+288 


1) GO TO 5 


O)VIPT(I.,1)=IPTCI, 1)4+2265 


O)IPT(1I,2)=IPTCI,2)4+2265 


-O)IPTCI,3)=IPTCI,3)+288 


O)IPT(CI.4)=IPT(I,4)+288 


O)IPT( I. 4} =IPTCI, 1)+288 
-O)IPT(I.3)=IPT(1I,3)+288 
-O)VIPT(I.4)=IPT(1I,4)+288 


ORDER ON FILE, NOT BY GRID 


SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SLVH 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 


Oana 


REORDERS POINT NUMBERS FOR READ 


32.162, 47.35, 
ew tOd. -1OG. 108.4 41055 
AAS = 61505., 1955-6 1602. 165.5 


10 IR=0 
DO 13 K=1,NP 
po 13 L=1,4 
MP=IPT(K,L) 
IF(MP.LT.1) GO TO 13 
41 II=K 
JU=EL 
DO 12 I=1,NP 
DO 12 J=1,4 
IF (IPT(1,JU).LT.1) GO TO 12 
IFC IPT(I.JU).GT.3490) GO TO 12 
IF(IPT(I,U).GE.MP) GO TO 12 
I1=1J 
JJ=J 
MP=IPT(I,J) 
42 CONTINUE 
IFCIPT(II,JJ).GT.3490) GO TO 14 
IR=IR+1 
IREAD(IR,1)=I11 
IREAD(IR,2)=Ju 
IREAD(IR,3)=IPTCII,JJ) 
IPT(II, uJ) =IPT(II,JJ)4+9000 
MP=IPT(K,L) 
IF(MP.GT.3490) GO TO 13 
GO To 11 
13 CONTINUE 
14 RETURN 
END 
SUBROUTINE STOATM(Z,T,P,D) 
DIMENSION 28(49),1TMS(49),WMS(49),PS(49) 
DATA (ZS(1),1=1,49)/0., 11.019, 20.063, 
& 54.413,71.802,86. ,91.,94.,97.,100 
&® 112.,115.,120.,125.,1390.,135.,140 
& 170., 180., 190., 210., 230., 265. 


& 500. ,550. ,600. ,650. , 700. ,750. , 800 


, 300., 9350., 400., 450., 
.,850. ,900. ,950. , 1000. / 


DATA (TMS(1),1=1,49)/288.15,216.65,216 65,228 .65,270.65,270.65, 

& 214.65, 186.95, 186.87, 187.74,190.40, 195.08, 202.23,212.89,223.29, 
& 240.0,264.0,300.00,360.00,417.23,469.27,516.59,559.63,598.78, 
& 634.39, 666.80, 696.29, 723.13, 747.57, 790.07, 825.31. 
& 878.84, 915.78, 955.20, 976.01, 990.06, 995.83, 998.22, 
&® 999.24, 999.67, 999.85, 999.93, 999.97, 999.99, 999.99. 
&® 1000., 1000., 1000., 1000./ 


DATA (WMS(1L),1=1,.49)/28.9644,28 .9644,28.9644,28.9644,28 .9644, 


28 
27 
24 
22 
15 
& 6. 
DAT 
& .6 


Ce Se Pe 


9644 28.9644 ,28.9522,28 .889,28.783,28.620,28.395,28. 104, 
_765,27.521,27.268,27.020,26.680, 26.205,25.803, 25.436,25.087, 
.749, 24.422, 24.103, 23.792, 23.488, 23.192, 22.902, 

.342, 21.809, 20.825, 19.952, 18.688, 17.726, 16.7935, 

.984, 15.247, 14.330, 13.092, 11.505, 9.718, 7.998, 

579, 5.543, 4.849, 4.404, 4.122, 3.940/ 

A (PS(1I),1=1.49)/1019.25, 226.32, 54.7487, 8.68014, 1.10905, 
6938, .039564, 3.7338E-3, 1.5381&-3, 9.OS560E-4, 5.3571E-4, 


SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
SORT 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 


795 


76 


12 
13 


Bi 


31 


32 


33 


GO GT BT BD @ we we 


.2011E-4,1.9742E-4,1.2454E-4,9.3188E-5,7.1042E-5,5.5547E-5, 


.1908E-10, 2.2599E-10, 1.7036E-10, 1.3415E&-10, 1.0873E-10. 
.981GE-11, 7.5138E-11/ 

IF(Z2.LT.O.) GO TO 81 

RO=6356.766 

GO=S . 80665 

WMO=28 .9644 

RS=8314.32 

ZM=2* 1000. 

ROM=RO* 1000. 

IF(Z.GE.86.) GO TO 6&6 

DO 3 I=1,7 
IF(ZSC(I).LE.Z.AND.2.LT.ZS(1+1)) GO TO 5 
CONT INUE 

ZL=RO*ZS(1)/(RO+ZS(1)) 
ZU=RO*ZS(1+1)/(RO+ZS(1+1)) 

ZLM=ZL* 1000. 

ZUM=2U* 1000. 

WM=WMO 

HT=(RO*Z)/(RO+Z) 

HM=HT* 1000. 

G=(TMS(1+1)-TMS(1))/(ZU-ZL) 

GM=G* .001 

IF(G.LT.O..OR.G.GT.0O.) GO TO 12 
P=PS(I)*EXP(-(GO*WMO* (HM-ZLM))/€RS*TMS(I)))* 100. 
GO TO 13 


3 

4.0096E-5, 2.5382E-5, 1.7354E-5, 1.25054E-5, 9.3568E-6, 
7.2028E-6, 5.6691E-6, 4.5422E-6, 3.6990E-6, 9.0395E-6, 
2.5278E-6, 2.12tOE-6, 1.5271E-6, 1.1266E-G6, 6.4756E-7, 
3.9276E-7, 1.7874E-7, 8.7704E-8, 3.4498E-8, 1.4518E-8, 
6.4468E-9, 3.0236E-9, 1.5137E-9, 8.2130E-10, 4.8865E-10, 
3 

8 


P=PS(1)*((TMS(1)/(1TMS(1)+G*(HT-ZL)))**((GO*WMO)/(RS*GM)))*100. 


T=TMS(1)+G*(HT-ZL) 
GO TO 25 
DO 7 1=8,48 
IF(Z2S(1).LE.2.AND.2.LT.ZS(1+1)) GO TO 8 
CONT I NUE 
1=48 
IF(Z.LE.1000.)GO TO 8 
T=O. 
P=0. 
D=O. 
RETURN 
IF(I.NE.8)GOTO 31 
T = T™MS(9) 
GOTO 39 
IF(I.LT.16.0R.1.GE.19)GOTO 32 
T = 240.+12.0*(2-110.0) 
GO TO 39 
IF(I.GE.19)GOTO 33 
T = 263.1905 - 76.3232*SQORT(1. - ((2-91.)/19.9429)**2) 
GOTO 39 
XI = (2-120.)*(RO + 120.)/(RO + Z) 
T = 1000. - 640. *EXP(-O.01875*XT) 


39 


24 


J=zil 

IF(I.EQ.48)J = I - 1 
ZO = ZS(d) 

Z1 = ZS(uJU+t) 

Z2 = 2S5(J+2) 


WMA=WMS(U)*(2-Z1)*(2-Z2)/(€Z0-Z1)*(Z20-Z22)) + WMS(U+1)*(2Z-Z0) 
& *(2-22)/((Z1-2Z0)*(2Z1-22)) + WMS(U+2)*(Z2-20)*(2-Z1)/ 
& ((22-20)*(22-21)) 
ALPO ALOG(PS(J)) 
ALP 4 ALOG(PS(J+1)) 
ALP2 = ALOG(PS(JU+t2)) 
ALPA = ALPO*(7-2Z1)*(2-22)/(€(Z0-21)*(Z0-2Z2)) + ALP1*(Z-20) 
& *(2-22)/((Z1-Z0)*(21-22)) + ALP2*(Z-20)*(Z-21)/ 
& ((22-Z0)*(2Z2-2Z1)) 
ALPB = ALPA 
WMB = WMA 
IF(I.€Q.8.0R.1.EQ.48)GOTO 24 
J=zrd- 1 
ZO = 2S(uU) 
Z1 = ZS(uJU+t) 
22 = 2S(J+t2) 
ALPO ALOG(PS(J)) 
ALP 1 ALOG(PS(uJ+tt)) 
ALP2 ALOG(PS(UJU+2)) 
ALPB = ALPO*(2-21)+(2-22)/((Z0-Z1)*€Z0-Z2)) + ALPIt*(Z-Z0) 
& *(2-22)/(€(2Z1-720)*(21-Z2)) + ALP2*(2-270)*(2-21)/ 
& ((Z22-20)*(2Z2-Z1)) 
WMB=WMS(JU)*(2-Z1)*(Z-22)/((Z0-21)*(Z0-22)) + WMS(J+1)*(Z-Z20) 
& *(2Z-22)/((Z1-Z0)*(21-22)) + WMS(JU+2)*(Z-20)*€2-Z1)/ 
&® ((22-70)*(22-21)) 
P= 100.*EXP((ALPA+ALPB)/2. ) 
WM = (WMA+WMB)/2. 
25 D=(WMtP)/(RS#T) 
26 RETURN 
END 
FUNCTION TEMP ( ALT , TX , Ti, 7T3', T4 , AD ) 


lou 


STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STOA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDA 
STDBA 
STDA 
STDA 
TEMP 
TEMP 


CR ERE REET RTT TRE TE EEE RTE OR TR RE Rk TE MD 


Cc ** 
Cc** 
C** 
Cc** 
Cee 


‘TEMP’ CALCULATES THE TEMPERATURE AT ALFIVTUDE ALT USING EQUATION 
(10) FOR ALTITUDES BETWEEN 90 AND 125 KM AND EQUATION (13) 
FOR ALTITUDES GREATER THAN 125 KM , FROM SAO REPORT 313. 


**TEMP 
++ TEMP 
++ TEMP 
** TEMP 
*+TEMP 


CER RRR ETE ER ERT OT OR KEY KE OR KO RR ER Rr RR KO KR RK ET KR TEMP 


PARAMETER (BB = 4.5E-6) 


U = ALT - 125. 
IF ( U .GT . oO. ) THEN 

TEMP = TX4+A2*ATAN(T1*U*(1.48B*(Ut*%2.5))/A2) 
ELSE 


TEMP ™ + Tt * U + TD * (UY*3) + T4 * (U**4) 


Hl 


END IF 


TEMP 
TEMP 
TEMP 
TEMP 
TEMP 
TEMP 
TEMP 
TEMP 
TEMP 
TEMP 





77 


78 


END 
SUBROUTINE TINF(F10,F10B,GI,XLAT,SDA,SHA,DY,11,TE) 


TEMP 
TINE 
TINE 


CUR ER ERE REE EE EEE ETE TT ETT EEE ETHER EEE EEE Ee HEE KE ETT INE 


C** SUBROUTINE ‘TINF’ CALCULATES THF EXOSPHERIC TEMPERATURE ACCORDING 


ct* TOL 
C+** 


. JACCHIA SAO 313, 1970 


ce F10 = SOLAR RADIO NOISE FLUX ( X E-22 WATTS / M2 ) 
C** FIOB= 162-DAY AVERAGE F10 


Ce GI 
Cz XL 


= GEOMAGNETIC ACTIVITY INDEX 
AT= GEOGRAPHIC LATITUDE AT PERIGEE ( IN RAD ) 


cr SDA = SOLAR DECLINATION ANGLE ( IN RAD ) 
c* SHA = SOLAR HOUR ANGLE 


i DY 
Cet 11 
C** ~ RE 
crt 
C.F 
Ct 
cf 
Ctr 
cr 


- 0 / Y ( DAY NUMBER / TROPICAL YEAR ) ; 1 
= GEOMAGNETIC EQUATION INDEX ( 14--GI=KP , 2--GI=AP ) 
= DIURNAL FACTOR KP, F10B, AVG 


CONSTANTS -- C = SOLAR ACTIVITY VARIATION 
-- BETA , ETC = DIURNAL VARIATION 
-- D = GEOMAGNETIC VARIATION 
-- £ = SEMIANNUAL VARIATION 


** TINE 
** TINE 
** TINE 
** TINE 
** TINE 
**T INE 
¥* TINE 
** TINE 
** TINE 
** TINE 
** TINE 
Sa LNE 
eT TNF 
** TINE 
** TINE 
** TINE 
** TINE 
** TINE 


CURT EEE RE ERR RT ET ET TERT REEDS RR RTT HEH eH TING 


Cc 
c** SOLA 
C 


C 


PARAMETER (PI = 3.14159265 . TP] = 6.28318531) 
PARAMETER (XM = 2.5 , XNN = 3.0) 


ARE SOLAR ACTIVITY VARIATION VARIABLES 

PARAMETER (C1 = 383.0 , C2 = 3.32 , C3 = 1.80) 

ARF GEOMAGNETIC VARIATION VARIABLES 

PARAMETER (D1=28.0,02=0.03,D3=1.0,D0D4=100.0,D5 = -0.08) 
ARE SEMIANNUAL VARIATION VARIABLES 

PARAMETER (£1=2.41,6&2=0.349,E3=0.206,£4=6.2831853) 
PARAMETER (E5=3.9531708,E6=12.5663706,E7=4.3214352) 
PARAMETER (E8 = 0.1145 . £E9 = 0.5 , E10 = 6.2831853) 
PARAMETER (E114 = 5.9742620 , E12 = 2.16) 


PARAMETER (BETA=-0.6457718,GAMMA=O0. 7504916 ,P=0. 1047198) 
PARAMETER (RE = 0.31) 


R ACTIVITY VARIATION 


TT = Ct + C2 * FIOB + C3 * ( FIO - FIOB ) 


C** DIURNAL VARIATION 


GC 


ETA = ©.5 * ABS ( XLAT - SDA ) 
THETA = 0.5 * ABS ( XLAT + SDA ) 
TAU - SHA + BETA + P * SIN ( SHA + GAMMA ) 
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TINF 
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TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINE 
TINF 
TINF 
TINF 
TINF 


OOnNNA UO AWN =O 


IF ( TAU. GT. PT ) TAU = TAU - TPT 
IF ( TAU. LT.-PI ) TAU = TAU + TPI 
A1 = ( SIN ( THETA ) )}**xXM 
A2 = ( COS ( ETA ) )**xXM 
A3 = ( cos ( TAU / 2. ) )**xXNN 
B1 = 1.0+ RE * Af 
B2 = ( A2 - Ai) / Bt 
TV = Bi * ( ¢. + RE * B2 * AD) 
TL = TC * TV 
C 
C** GEOMAGNETIC VARIATION 
Cc 
IF ( 11.€0.1 ) THEN 
TG = D1 * GI + D2 * EXP(GI) 
ELSE 
TG = D3 * GI + 04 * ( 1 
END IF 
C 
C** SEMIANNUAL VARIATION 
Cc 
G3 = 0.5 * ( 1.0 + SIN ( E10 * DY + E11 
G3 = G3 ** E12 
TAU = DY + E8 * ( G3 - EQ ) 
Gi = E2 + £3 * ( SIN ( E4 * TAU + ES ) ) 
G2 = SIN ( EG * TAUI+ E7 ) 
TS = Ei + FIOB * Gi * G2 
C 
C** EXOSPHERIC TEMPERATURE 
Cc 
TE = TL + TG + TS 
RETURN 
END 
SUBROUTINE TME ( MN , IDA , IYR , IHR , 
SDA , SHA , DD, DY ) 


MIN 


- EXP ( DS * GI ) ) 


) ) 


» XLAT , XLNG , 


TINF 
TINE 
TINF 
TINF 
TINE 
TINF 
TINE 
TINE 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINF 
TINE 
TINE 
TINE 
TINF 
TINF 
TINE 
TINE 
TINF 
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TINF 
TINF 
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TINF 
TINF 
TINF 
TINF 
TME 

TME 

TME 
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-180,+180) 


C** ‘“TME’ PERFORMS THE CALCULATIONS OF THE SOLAR DECLINATION 
C** ANGLE AND SOLAR HOUR ANGLE. 

C** 

C** INPUTS: MN = MONTH 

C** IDA = DAY 

C*+* IYR = YEAR 

C** IHR = HOUR 

c** MIN = MINUTE 

C** XLAT= LATITUDE ( INPUT-GEOCENTRIC LATITUDE ) 
C** XLNG= LONGITUDE (INPUT-GEOCENTRIC LONGITUDE, 
C** 

C** OUTPUTS: SDA = SOLAR DECLINATION ANGLE (RAD) 

C** SHA = SOLAR HOUR ANGLE (RAD) 

C** DD = DAY NUMBER FROM 1 JAN. 

C+ DY = OD / TROPICAL YEAR 


**TME 
**TME 
**TME 
**TME 
**TME 
** TME 
**TME 
**TME 
**TME 
**TME 
**TME 
**TME 
**TME 
**TME 
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20 


PARAMETER (YEAR = 365.2422) 
PARAMETER (A1 99.6909833 , A2 
PARAMETER (A3 0.00038708 , A4 O.250684477) 

PARAMETER (B1 0.0172028 . B2 0.0335 , BI = 1.407) 
PARAMETER (PI 3.14159265 . TPI = 6.28318531) 

PARAMETER (PI2 = 1.57079633 , PI32 = 4.71238898) 

PARAMETER (PI180 = 0.017453293) 


36000. 76892) 


DIMENSION IDAY(12) 
DATA IDAY / 31,28,31 ,30,31,30 ,31,31,30 .31,30,31 / 


XLAT = XLAT / 57.29577951 
YR = IYR 


IF ( MOOCITYR,4) .EQ. O ) THEN 
IF ( MOD(TYR,100) .NE. O ) IDAY(2) = 29 
ELSE 
IOAY(2) = 28 
END IF 
ID = 0 
IF ( MN. GT. 4 ) THEN 
DO 20 I = 41 +, MN-14 
ID = ID + IDAY(I) 


CONTINUE 
FND IF 
ID = ID + IDA 
DD = ID 
DY = DD/YEAR 


COMPUTE MEAN JULIAN DATE 


XMUD = 2415020. + 365. * ( YR - 1900. ) + OD 
+ FLOAT ( ( IYR - 1901 ) / 4 ) 


COMPUTE GREENWICH MEAN TIME IN MINUTES GMT 
XHR = IHR 
XMIN = MIN 
GMT = 6O * XHR + XMIN 
FMJD = XMJD - 2435839. + GMT / 1440. 


COMPUTE GREENWICH MEAN POSITION - GP ( IN RAD ) 


( XMUD - 2415020.5 ) / ( 36525.0 ) 
AMOD ( At + A2 * XU + AB * XJ * XU + AS * GMT , 3EO. ) 


XJ 
GP 


'] 


COMPUTE RIGHT ASCENSION POINT - RAP ( IN RAD ) 
1ST CONVERT GEOCENTRIC LONGITUDE TU DEG LONGITUDE - WEST , + EAST 


IF ( XLNG .GT. 180. ) XLNG = XLNG - 360. 
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RAP = AMOD ( GP + XLNG , 360. ) 


Cc 
C** COMPUTE CELESTIAL LONGITUDE - XLS ( IN RAD ) - - ZERO TO 2PI 
Cc 

Yi = B1 * FMUD 

Y2 = 0.017202 + ( FMJD - 3. ) 

XLS = AMOD ( Yi + B2 * SIN(Y2) - B3 , TPI ) 
C 


C** COMPUTE SOLAR DECLINATION ANGLE - SDA ( IN RAD ) 


B4 = PI18O * ( 23.4523 - 0.013 * XJ ) 
SDA = ASIN ( SIN ( XLS ) * SIN ( B4 ) ) 


C 
C** COMPUTE RIGHT ASCENSION OF SUN - RAS ( IN RAD ) - - ZERO 10 2PI 
C 
ARG = TAN ( SDA ) / TAN ( Ba ) 
IF ( ARG .GT. 1.0 ) ARG = 1.0 
IF ( ARG .LT. -1. ) ARG = -1.0 
RAS = ASIN ( ARG ) 
C 
C** PUT RAS IN SAME QUADRANT AS XLS 
C 
RAS = ABS ( RAS ) 
TEMP = ABS ( XLS ) 
IF ( TEMP.LE.PI .AND. TEMP.GT.PI2 ) THEN 
RAS = PI - RAS 
ELSE IF ( TEMP.LE.PI32 .AND. TEMP.GT_PI ) THEN 
RAS = PI + RAS 
ELSE IF € TEMP.GT.PI32 ) THEN 
RAS = TPI - RAS 
END IF 
IF (€ xtS. LT. O. ) RAS = -RAS 
C 
C** COMPUTE SOLAR HOUR ANGLE - SHA ( IN RAD ) - - 
Cc 


SHA = RAP * PI]180 - RAS 
IF ( SHA.GT.PI ) SHA = SHA - TPI 
If ( SHA.LT.-PI ) SHA = SHA + TPI 


RETURN 
END 
SUBROUTINE WIND 
COMMON/WINCOM/RHO,FCORY ,DX5,DYS,PX,PY,U,V.T,TX,TY, 
$ DU,DV,P,UPRE,VPRE ,DUPRE ,DVPRE 
COMMON/ IOTEMP/DUM1(6),PHI,O0UM2(11).MN,0M2A(5),G.R,H,PHIR, 
& THETR,DUM3(15),FLAT,DUMMY( 18) 
COMMON/CHIC/DUM(18), I WSYM,USH, VSH,OUS!H!,DVSH 
ABSPHI = ABS(PHIR) 
BETA = 1.458E-6 


TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
TME 
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on DD oO 
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20 


25 


31 


40 


9S 


SVAL = 110.4 
SHZ1=25. 
S$HZ2=120. 
IF (ABSPHI .LT.0.017453293*FLAT)GO TO 40 


IF(RHO.GT.O.0.AND.T.GT.0.0.AND.ABS(FCORY) .GT.O. )GO TO 20 


IF (RHO.GT.O..AND.ABSPHI .GT.O.) GO TO 20 
U= O. 
V =O. 
DU=O. 
DV=0. 
IF(ABS(FCORY).LE.O.)GO TO 31 
RETURN 
FCORX = FCORY*OX5/DYS 
IF(H.GE.90.)GOTO 25 


U = - PY/(FCORY*RHO) 

V = Px /( FCORX*RHO) 
DU = -(G*TY)/(FCORY*T) 
DV = (G*TX)/(FCORX*T) 
Go To 31 


VISC = BETA*(T**1.5)/(T+SVAL) 
VLS = 5.3 + 0.0622*(H**1.5) 
IF(VLS.GT.H)VLS = H 
VISCFAC = VISC/(1t.OEG*RHO'TVLS**2) 
CORIOL = FCORY/DY5S 
DPDX = PX/(DX5*RHO) 
DPDY = PY/(DYS5*RHO) 
DENOM = CORIOL**2 + VISCFAC**t2 
U = (-CORIOL*OPDY - VISCFAC*DPDX)/DENOM 
Vv = (CORIOL*DPDX - VISCFAC*DPDY )/OENOM 
DTDX = G*TX/(DxX5*T) 
DTIDY = G*TY/(DY5*T) 
DU = (-CORIOL*DTDY - VISCFAC*OTDX)/DENOM 
DV = (CORIOL*DTDX - VISCFAC*DTDY )/DENOM 
IF(H.GT.SHZ4.AND.H.LT.SHZ2)GOTO 39 
IF (ABSPHI .GE.0.017453293*FLAT )RETURN 
CONT INUE 
U=UPRE 
V=VPRE 
DU=DUPRE 
DV=DVPRE 
IF(H.GT.SHZ1.AND.H.LT.SHZ2)GO TO 99 
RETURN 
IH=INT(H) 
IFCIH.LT.SHZ1+5) GOTO 130 
IF(IH.GE.SHZ2-30)GOTO 140 
U=USH 
V=VSH 
DU=DUSH 
OV=DVSH 
RETURN 


C...LOW ALTITUDE FAIRING 
130 


FACS=(H-SHZ1)/5. 
FACG=1.-FACS 
U=FACG*U+FACS*USH 
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V=FACG*V+FACS*VSH 
DU=F ACG*DU+FACS*DUSH 
DV=F ACG*OV+FACS*DVSH 
RETURN 


C...HIGH ALTITUDE FAIRING 


140 


FACS=(H+30.-SHZ2)/30. 


FACG=1.-FACS 
U=FACS*U+FACG*TUSH 
V=FACS*V+FACG*VSH 
DU=FACS *0U4+FACG*DUSH 
DV=FACS*OV+FACG*tDVSH 
RETURN 

END 


WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
WIND 
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2. JCL CODE 


GRAM-90 on MSFC’'s IBM 3090 computer. 
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//CWEJ412T JOB (665551500466), ‘JEFFRIES/5449136'. 

// MSGLEVEL=1,MSGCLASS=X,CLASS=A, REGION=2048K 

/*JOBPARM F=S9700.R=412 

//GRM9OTST EXEC VSF2CLG.FVLNSPC= ‘3200, (200,50), 

// PARM.FORT=’SDUMP( ISN) ’ 

//? a Nt) ss tress Smet Clay as sh cs mere i Sa a ae a eam ne eo ea ee 
//* THIS RUN STREAM WILL COMPILE THE VS FORTRAN SOURCE PROGRAM 

//* WHICH IS MEMBER GRAM9O IN THE CATALOGED PARTITIONED DATA SET 
//* CWEJ412.ATMOSMDL FORT - A VERSION OF THE GRAM. 


ff 

//* LAST MOD: 10 SEP 1990 

//* bwawtoLtwwtiuweae ac Loee Skee oe Se Meee See = See ee Se eee SSeS See ee 
//FORT.SYSIN DD DSNAME =CWEJU412.ATMOSMDL . FORT (GRAMSO) ,OISP=SHR 
18 ied 

//* INPUT UNIT FOR DATA LINES(SEE INPUT DATA BELOW). 

//* 

//GO.FTOSFOO1 DD DONAME=SYSIN 

//* 

/f/* OUTPUT PRINT UNIT. 

fi 

//GO.FTOGFOO! DD SYSOUT=x 

J f% 


//* OUTPUT DATA SET. IF USED, IT WILL STORE THE RAW DUTPUT DATA 
//* IN THE CHOSEN DATA SET FOR FUTURE PLOTTING/OTHER HANDLING. 

i fad CHANGE AS REQUIRED. 

i es 

//*O0.FTO7TFOO1 DD DSN=CWEU412. PLOT .DATA(GRMSOTST) ,DISP=SHR 
//*0.FTOTFOO1 DD DSN=CWEJU412.STSMET .DATA(G8STS28A ) ,.DISP=SHR 

//* 

//* INPUT DATA SETS(1 SCIDAT AND 1 4-D SETS): 

//* 

//* DSN CWEJU412.GRAMODSO.DATA(SCIDATS) INPUTS SCIDAT DATA. 

Lf 

//GO.FTO3FOO1 DD DSN=CWEJ412.GRAMOD9O.DATA(SCIDATS) ,DISP=SHR 

ifs 

//* DSN CWEU412.GRAMODSO .DAT4DXX, A SERIES OF DIRECT ACCESS (DA) FILES, 
ie bd INPUTS 4D DATA, WHERE XX IS THE NUMBER OF THE DESIRED MONTH. 
ii CHANGE THE DSN NAME XX AS DESTRED (OR IN ACCORDANCE WITH AN 
1 fs ALTERNATE DA FILE NAME): 

//* JAN. FEB. MAR, APR.MAY, JUN, JUL. AUG.SEP.OCT.NOV.DEC. 


i /* O1, 02, 03, 04, O05, O06. 07, 08. O09. 10, 11, 12, 

Fe de BUT THE FOLLOWING "GO" NAME MUST BE COMPATIBLE WITH THE "OPEN" 
//* STATEMENT FILE NAME IN THE FORTRAN PROGRAM. 

PS 


//GO.DAT4DOS ODD DSN=CWEJ412.GRAMODSC.DAT4DOS, 

//  DCB=(RECFM=F ,DSORG=DA, LRECL= 1269) .DISP=SHR 

1 es 

//* THE FOLLOWING UNITS ARE ASSIGNED TO INTERNAL SCRATCH FILES. 
//* LRECL OF 852 = 213 (SIZE, 1.£. NUMBER OF WORDS, ALLOCATED FOR INPUT 
//* FILE IN SUBROUTINE GRID4D) X 4 (BYTES/WORD), AND BLKSIZE = LRECL 
f/* XAG; 

Std 

//GO.FT12FOO1 DD UNIT=SYSDA,SPACE=(TRK,(20,.20)). 

J / DCB=(RECFM=FB,LRECL=852 ,BLKSIZE=8520) 


The following is a sample listing of the Job Control Language (JCL) code needed to run the 


//GO.FT13FOO1 OD UNIT=SYSDA,SPACE=(TRK,(20,20)), 
// OCB=(RECFM=FB,LRECL=852,BLKS1ZE=8520) 


ie ia 

//* INPUT DATA: 

//* LINE 1: INITIAL ALT(KM),INITIAL LAT,INITIAL WEST LON, 

//* F10.7,MEAN F10.7,AP,MONTH,DAY,YR(2 DIGITS), 

//* HR(24 HR CLOCK) ,MIN,SEC,.LAT INCREMENT,.W. LON 

//* INCREMENT ,ALT INCREMENT ,MAX #4 POSITIONS, TIME 

//* INCREMENT , TRAJECTORY OPTION,OUTPUT OPTION( SHOULD 

//* EQUAL 7 IF ACTIVATED), MIN GEOSTROPHIC LAT. 

i f* NOTE: FOR STS SUPPORT: 

f/* INIT ALT=122.8 KM (403K FT) 

//* ALT INCR=.9 KM (3K FT) 

//* MAX # POS=138 (TO SFC) 

//* USUALLY IN AN OPERATIONAL SITUATION, THE ONLY 

//* DATA ITEMS WHICH NEED TO BE CHANGED ARE THE: 

//* F10.7,MEAN F10.7,AP,MONTH,DAY,YR,HR, AND MIN. 
//* LINE 2: INPUT UNITS FOR DATA FROM SCIDAT (NMC,.GROVES/STA PERT. 
//* RANDOM ,P/PW/CS/CL,QB0,SP),4D,RANOOM OPTION, QBO 

{//* OPTION,FIRST RANDOM #,4D/P/D/T SCRATCH UNIT, 

//* NMC SCRATCH UNIT, AND DIAGNOSTIC OUTPUT SWITCH, IPRT 
//* (O ACTIVATES OUTPUT, 1 DE-ACTIVATES IT). ALSO 

//* NLIMIT(MAX NUMBER OF TEST VIOLATIONS IN 4D DATA - 12 
//* TO 25). 

/f/* LINE 3: INITIAL PL,PS,OL,0S,TL,TS,UL,US,VL,VS,RANDOM 

Ee Migs PERT SCALING FACTOR--INCLUDE ONLY IF RANDOM 

/{* OPTION=1(IF 2, RANDOM PERT NOT COMPUTED - SAME 

1 /* FOR QBO OPTION). 

//* LINES 4 OR MORE: ONLY REQUIRED IF TRAJECTORY OPTION .NE. 

//* ZERO. EACH LINE HAS TIME(SEC),ALT(KM).LAT,W.LON. 

a LAST LINE HAS ANY NEGATIVE ALTITUDE VALUE. 

Ps * 


//GO.SYSIN DO * 

600.0,28.45,80.53, 150.,150.,15.,.9,1,89,00.00,0, 0., 0.,20.,31,0,0,0, 18., 
3, 4, Ty, Ve 1.92. 139.4. 12, 

O:,0.70..:0.,07.0-.0:102.0:.0:%..41 0; 

/* 

// 
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3. OUTPUT 


The following is the sample output, from MSFC’s IBM 3090 computer, resulting from the 
input submitted through the foregoing JCL listing. As in a similar example listed in Volume | 
(Justus et al. 1990), the first 18 lines merely repeat the input data for ease in comparison and 
reference. 
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4. 4-D CONVERSION 
4.1 Description 


The world-wide meteorological data set developed for the 4-D model by Allied Research 
Associates (Spiegler and Fowler, 1972) was originally binary tapes labeled WWIA-WW3A. The 
files on these binary tapes were split up by many users into individual monthly binary or ASCII 
sequential files stored on their computers. These files, one per month, have now been converted to 
ASCII files which can be read by random access (or “direct access”) I/O operations in the new 
GRAM-90. 


Within each file are 3,490 records representing the values at individual grid points of latt- 
tude and longitude. These points are grouped into three grid areas: 288 points on the Northern 
Hemisphere equatorial (EQN) grid: 1,977 points on the Northern Hemisphere (National Meteorolog- 
ical Center, NMC) grid; and 1,225 points on the Southern Hemisphere (SH) grid. On the NMC 
grid. the data were computed at NMC points and are stored in the order given by the NMC grid 
table given in the SCIDAT9 data file which concludes this volume. For the other two grid areas. 
the data are given at 5° Jatitude-longitude intersections westward from the Greenwich Meridian to 
5° E. The EQN grid covers the latitudes from 0° to 15° N. with points occurring in the following 
order: 1-4 =Ilong. 0. lat. 0. 5, 10. 15: 5-8 =long. 5° W., lat. 0, 5. 10, 15: ...285-288 = long. 3° 
E.. lat. 0. 5S. 10. 15. The SH grid contains all data from 5° S. to the South Pole as follows: 
| =South Pole, 2-8 = long. 0. lat.—5 to -85: 19-35 =long. 5° W.. lat. -5 to -85. ...1,209-1.225 
=long. 5° E.. lat. —S to —85. It should be noted that the South Pole is given only once. as the 
first point of the SH data set. 


Each record consists of 1,268 ASCII characters (formatted 208F6.2, 2F4.2. 3F4.0). The 
first 208 parameter values contain the thermodynamic data for a latitude-longitude grid point. The 
last five variable values are identifiers. The data are arranged by altitude level for each parameter: 
thus. the first 26 values contain the pressure means from the surface to 25 km and the next 26 
values contain the pressure variances for the same levels. This pattern continues for the 26 levels 
of temperature means and variances, moisture means and variances, and density means and vari- 
ances. Temperatures are in degrees Kelvin (K), pressures are in millibars (mb) (multiplied by 100 
for use in GRAM-9O, to give Newtons per square meter (N-m--)), and densities are in grams per 
cubic meter (g-m™*) (divided by 1,000 for use in GRAM-90, to give kilograms per cubic meter 
(kg-m“)). GRAM-90 does not use the moisture means and variances. 


Variables 209 and 210 contain the latitude and longitude of the data grid point (in 1Oths of 
a degree). The latitude is always positive (since the Southern Hemisphere is identified by grid), 
and the longitude 1s always west. 


The last three values (variables 211-213) contain three identifier parameters. The first of 
these is the homogeneous moisture region in which the point lies, the second is the point number, 
and the third is the month. It should be noted that the 4-D points are numbered within the grid that 
contains them, and not by their location in the file. Thus, the point numbers run from 1-288 (EQN 
grid), 1-1.977 (NMC grid), and 1-1,225 (SH grid), not from 1-3,490. 
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The following JCL and FORTRAN code sequence is one example which has been used on 
MSFC’s IBM 3090 computer to reconfigure the 4-D sequential data sets needed by previous 
GRAM versions into the direct access data sets required by this version of the GRAM-90. In this 
example, conversion is made from a sequential ASCII format version of a single-month data set to 
the ASCII direct access format. It should be noted that due to several previously identified errors 
in the April, June, and October 4-D data bases, most users may have incorrect versions of those 
months, which must be corrected or replaced before conversion to enable correct functioning. 
Corrected copies of or corrections to the data bases may be obtained from Mr 
ES44, NASA MSFC, Huntsville, AL 35812. 


4.2 Code 
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//CWEU412C JOB (6ES551500466). ‘JEFFRIES/5449 136’ 

// MSGLEVEL=1.MSGCLASS=X .CLASS=A 

/*JOBPARM £=9700,R=412 

//***THE FOLLOWING STEP PREALLOCATES A DIRECT ACCESS DATA SET WHICH IS 
//***LATER INVOKED, LINKED TO A UNIT AND MANIPULATED. 

//***LATEST MOD: 22 MAY 1990 

//STEPO EXEC PGM=IEFBR14 


//O0D1 OD OSN=CWEJU412.GRAMODSO .NDAT4ADIO, 

//***NOTE: OSN MUST SBE CHANGED FOR EACH NEWLY ALLOCATED OATA SET, 

{fr MATCHED WITH THE ENTRY FOR FTOSFOO: BELOW, ANDO BE COMPATABLE 
ff/**? WITH THE SOURCE IN FTO7TFOO1 

Yi 4 DISP=(NEW,CATLG) .UNIT=SYSOA. 

ff SPACE=(1269.(3490.0)). 

//*"*NOTE: SPACE = BLKSIZE (].E€. # BYTES). » BLKS, # AOONL BLKS. 

ff DCB=(RECFM=F ,LRECL=1269 .O0SORG=DA) 


//***THE FOLLOWING STEP COMPILES. LINKS AND RUNS THE PROGRAM CONV4D 
//***WHICH IS INCLUDED IN THIS UCL STREAM. 
//CONV4D EXEC VSF2CLG,PARM FORT=’SDUMP(ISN)’ 
//FORT.SYSIN DO * 
PROGRAM CONV4D 


Cc READS IN 4-0 DATA TYPE IN MULTIRECORD FORMAT AND CONVERTS IT 
C TO A FORMAT SUITABLE FOR DIRECT (RANDOM) ACCESS. 
C LAST MOD: 22 MAY 1990 
C 
INTEGER DATA(213) 
C 


OPEN(8.ACCESS='DIRECT’ , FORM= ‘FORMATTED’ .RECL=1269) 
DO 10 [21,3490 
READ(7.8000.END=11) DATA 
8000 FORMAT(318.21(1018)) 
WRITE(8.BO10.REC=I1) OATA 
8010 FORMAT( 20BI6.514) 
10 CONTINUE 
11 CONTINUE 
END 
/* 
//GO.FTO7FOO!1 OO OSN=CWEU412.ATMOSMOL NDOATWW10( OCT). DISP=SHR 
//***THE FOREGOING IS THE SOURCE FILE OF THE DATA TO BE CONVERTED. 
//GO.FTO8FOO! OD OSN=CWEU412. GRAMODSO NDAT4010,O01SP2SHR 
eee FOREGOING IS THE TARGET FILE FOR THE CONVERTED DATA. 


. Dale Johnson, 


5. 4-D RECORD 


The following example depicts the sequence and form (but not the GRAM-90 format) of the 
values of the variables found in each individual data record of a monthly 4-D data set. The 
depicted data record is record number 3,348 of the April data set, a record in the Southern Hemi- 
sphere (SH) grid with NMC point number of 1,083, as previously described. 


99340 87590 77060 

67630 59180 51620 44890 38890 33560 28790 24710 21200 18140 
15550 13930 11430 9800 8400 7200 6170 5280 4530 3880 
3320 0 OQ 9011 7299 7672 BAQ7 9066 9344 9544 
9494 9171 809 1 5716 4032 2463 1905 649 293 113 
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6. SCIDAT9 DATA BASE 
6.1 Description 


NMC Grid Data. This data segment gives the 4-D Northern Hemisphere point number and 
the dual index for the corresponding NMC location. The NMC grid locations form an octagonal 
array, centered on the North Pole. The points are at square grid locations on the polar projection 
used for the NMC grid. A conversion between the latitude and longitude (treated as polar coor- 
dinates on the flat NMC grid plane) and the NMC grid indices (treated as Cartesian coordinates on 
the projection plane) is accomplished by a polar to Cartesian coordinate transformation, via 
equations programmed into the 4-D model. The NMC grid data in the SCIDAT9 file merely estab- 
lishes the equivalence between the sequential 4-D NMC point number and the 2-D x-y NMC grid 
point location. The NMC grid data constitute the first group of data in the SCIDAT9 file. The 
NMC grid data file contains 396 FORTRAN readable records with code “N” and I5 integers (A2, 
1517 format) in each record. 


Zonal Mean Data. The zonal mean data consist of 12 monthly sets of zonal mean values 
for pressure, density, temperature, and zonal wind, tabulated at 10° latitude intervals (from —90 to 
+90) and 5-km height increments (from 20 to 120 km) for each month. Prefix codes ZP. ZD. ZT, 
and ZU indicate pressure, density, temperature, and zonal wind, respectively. Each record contains 
the code, the month, the height in km, and the -90, -80, ...80, 90° latitude values of the parame- 
ter expressed as a four-digit integer, with an exponent common to all of the values in the field 
appearing at the end of the record. Thus, a value of 2,761 with an exponent at the end of the 
record of —6, would be the same as 2.761 x 10°° = 2.761 x 107%. Pressure data are in units of 
Nim-, density values are in kg/m". temperatures are in K, and zonal winds are in m/s. The zone 
mean data group contains 1.008 FORTRAN readable records, with the code and 22 integer values 
in each record (format A2, 14, 15, 1916, 14). 


Stationary Perturbations. The stationary perturbations are latitude-longitude dependent 
relative perturbations to be applied to the zonal mean values. Data for each of 12 months are given 
for both Northern and Southern Hemisphere latitudes. Prefix codes SP, SD, ST, SU, and SV indi- 
cate stationary perturbation values for pressure, density, temperature, zonal (eastward), or meridi- 
onal (northward) wind components, respectively. Each record contains the code, the month, the 
height in km, the latitude (-80 to +80) in degrees, and then 18 values of stationary perturbations, 
in per mil (%/10) tor the thermodynamic variables, and 0.1 m/s for the winds at longitude 180, 
160° W., 140° W., ...140° E., 160° E. The monthly mean value y,, for parameter y¥ at any latitude 
and longitude can be computed trom the zonal mean value z,, at the latitude, and the stationary 
perturbation s, (in per mil) at the latitude and longitude, by the relation 


Van = z.(1 ar Sy: /| 5000) 


Note that the stationary perturbation values at 90° latitude are always zero. The stationary perturba- 
tion data consist of 15,300 FORTRAN readable records, with a code and 21 integer values (format 
A2, 2175). 
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Random Perturbation Data. Random perturbation magnitudes (standard deviations) are 
latitude dependent only. Pretix codes RP, RD, RT, RU, and RV indicate random perturbation mag- 
nitudes in pressure, density, temperature, zonal wind, and meridional wind components, respec- 
tively. Each random perturbation record has the code, the month and height in km, followed by 19 
values of random perturbation magnitude, at 10° latitude increments (—90 to +90), followed by a 
common exponent value. These data give the relative standard deviations o, /p, Op /p, and of /T 
(in percent) for use in the random perturbation model. The code RU and RV data are similar, 
except that the wind perturbations are absolute deviations in m/s, and cover the height range 
0-200 km, whereas the RP, RD, and RT data cover the 20—-200-km range. Random perturbation 
magnitudes for the O-25-km altitude range are provided by the 4-D data base for the thermo- 
dynamic variables. The random perturbation data consist of 1,596 FORTRAN readable records with 
a code and 22 integer values (A2, 14, [5, 1916, [4) in each record. 


Large Scale Fraction Data. From daily difference analysis described in section 2 of 
Justus et al. (1980), the fraction of the total variance (o~ from the random perturbation data) con- 
tained in the large scale perturbations has been determined as a fraction of height and latitude. The 
SCIDATS file contains the annual average fraction (expressed as per mil) of total variance con- 
tained in the large scale. Large scale and small scale magnitudes, a, and o, are computed from the 
fractional data f, in per mil (code P for pressure, density, and temperature, or code PW for winds) 
by the relations 


a, = (Vf; 11,000) o, 


go, = (V1-f; /1,000) o, . 


where oy; is the total perturbation magnitude. The code P and code PW data groups each contain 
25 FORTRAN readable records, with code word P or PW, followed by 17 integer values (A2, 
| 717) on each record. 


Density-Velocity Correlations. Daily difference analysis (Justus et al., 1980) was also 
used to evaluate the cross correlations R,,, and R,,, for use in the velocity perturbation model, 
described in section 2 of part I. Both large scale and small scale values of the density-velocity 
correlations were evaluated, and are given in the SCIDAT9 data file (codes CL and CS) in per mil 
(1.e., divide by 1,000 to get correlations in the range —1 to +1). The code CS and CL data consist 
of 50 FORTRAN readable records, with code word CS or CL, followed by 12 integer values (A2, 


1217) in each record. 


Quasi-Biennial Oscillation (QBO) Data. The QBO data consist of height and latitude 
dependent amplitudes and phases for quasi-biennial variations in pressure (QP), density (QD), 
temperature (QT), and eastward and northward wind components (QU and QV, respectively). The 
amplitudes of the QBO wind components are in decimeters per second (0.1 m/s). The phases of 
all of the QBO parameters are measured in days after January 0, 1966, for the occurrence of the 
first maximum value. Since the period of the QBO variations is taken to be 870 days, the phases 
can vary from 0 to 870. Each QBO data record contains the code, the height in km, the amplitude 
and phase for 10° latitude, the amplitude and phase for 30° latitude, etc. out to the amplitude and 
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phase for 90° latitude. There are 80 FORTRAN readable records in the QBO data group. Each 
record contains 11 integer values, following a code word (QP, QD, QT, QU, or QV) in format 
A2, 1517. 

The following 1s a listing of the entire SCIDAT9 data base, which is further described in 
section 4.3 of volume I (Justus et al., 1990). 


6.2 DATA 
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